Immune Cytolytic Activity Stratifies Molecular Subsets of Human Pancreatic Cancer

Pancreatic adenocarcinoma (PDA) is the third most common cause of death from cancer, with an overall 5-year survival rate of less than 5%, and is predicted to become the second leading cause of cancer mortality in the United States by 2030 (1, 2). The American Cancer Society predicts that for the first time, more patients will die annually of pancreatic cancer than breast cancer beginning in 2016 in the United States (3). Recently approved combination chemotherapies for metastatic PDA modestly impact patient outcomes and durable remissions are rare (4, 5). Several recent studies have identified distinct genetic and transcriptional PDA tumor and stromal subtypes, which may present opportunities to identify individual patients likely to respond to targeted therapies (6–10). Immune modulation is a particularly attractive approach to treatment because of its potential to generate durable clinical responses in the proper setting (11, 12). Although single-agent immunotherapies targeting the immune checkpoint pathways PD1/PD-L1 and CTLA4 have shown striking efficacy in multiple tumor types, such approaches have failed to show clinical benefit in the overwhelming majority of patients with PDA (13, 14). Immunologically, PDA is characterized by a highly suppressive tumor microenvironment and a dense desmoplastic stroma (12, 15), and for most patients, there is scant intratumoral infiltration of effector T cells (6, 16). A small fraction of human PDA tumors do exhibit an immunogenic profile (6, 17), and there is provocative evidence that survival is improved in resectable PDA patients whose tumors have higher-than-average or unusual tumor T-cell infiltration (17–19). At present, the determinants of immune activation in PDA are poorly understood, providing little therapeutic guidance.

Data for tumor types available from The Cancer Genome Atlas (TCGA) were accessed in December 2015 and represent only untreated primary tumors (defined by the TCGA pathologist as “pancreatic ductal adenocarcinoma”; Supplementary Table S1). Each tumor sample is paired with a normal tissue sample providing a germline reference. The following tumor types (project code and n = sample size) were selected: kidney renal clear cell carcinoma (KIRC, n = 606), lung adenocarcinoma (LUAD, n = 116), cervical squamous cell carcinoma and endocervical adenocarcinoma (CESC, n = 309), lung squamous cell carcinoma (LUSC, n = 553), pancreatic adenocarcinoma (PAAD, n = 134), stomach adenocarcinoma (STAD, n = 418), head and neck squamous cell carcinoma (HNSC, n = 566), colon adenocarcinoma (COAD, n = 328), skin cutaneous melanoma (SKCM, n = 105 primary tumors), bladder urothelial carcinoma (BLCA, n = 427), esophageal carcinoma (ESCA, n = 196), liver hepatocellular carcinoma (LIHC, n = 371), thyroid carcinoma (THCA, n = 572), ovarian serous cystadenocarcinoma (OV, n = 309), glioblastoma multiforme (GBM, n = 169), and prostate adenocarcinoma (PRAD, n = 555). Manually curated DNA variant mutational annotation format and tumor mRNA expression were obtained for each disease type from Broad Firehose (http://gdac.broadinstitute.org). Raw DNA reads (.bam format) used for BMR calculation coverage estimates and HLA typing were accessed via the NCI Cancer Genomics Hub (https://cghub.ucsc.edu). GISTIC2.0 (20) individual copy number data for the PDA dataset was obtained from the TCGA Data Portal (https://tcga-data.nci.nih.gov).

Data were normalized following the method of Rooney and colleagues (21). Briefly, total raw read counts per gene were divided by the gene's maximum transcript length to represent a coverage depth estimate. Coverage estimates were then scaled to sum to a total depth of 1e6 per sample and can be interpreted as Transcripts Per Million (TPM; ref. 21). Gene set variation analysis (GSVA) was performed using the R/Bioconductor package “GSVA” because it implements a nonparametric unsupervised method to measure gene set enrichment across a dataset. The sample-wise enrichment score for a given gene set is calculated using a Kolmogorov–Smirnov (KS)-like random walk statistic. Statistical ranking for GSVA scores for the cytolytic index by the top decile and bottom quartile were defined as cytolytic-high and cytolytic-low, respectively. Unsupervised hierarchical clustering, using complete linkage with the distance metric equal to 1 minus the Pearson correlation coefficient, was also performed using the GSVA scores for each dataset. Differential gene expression analysis between cytolytic-high (top 10th decile cytolytic index) and cytolytic-low (bottom 25th quartile cytolytic index) across TCGA datasets was calculated using gene-level raw counts with the R/Bioconductor package “limma” with voom transformation with quantile normalization (22). Lowly expressed genes with less than 1 count per million in fewer than half of the samples in a dataset were excluded for differential gene expression analysis. Genes with BH-adjusted P values ≤ 0.1 were considered differentially expressed. All plots and graphs were generated using the R package “ggplot2” (23).

Significantly mutated genes (SMG) in cytolytic subtypes of PDA were calculated using the Mutational Significance in Cancer (MuSiC Genome Suite; ref. 24). MuSiC identifies SMGs with a significantly higher mutation rate than the background mutation frequency (BMR) for a given gene calculated across the entire sample population. The threshold for significance was a false discovery rate of 0.1. Mutational spectra across cytolytic subtypes of PDA were determined as described previously (10, 21). Somatic copy number alterations in each TCGA PDA sample were counted by taking the sum of segment mean changes ≥ 0.6 and ≤ −0.4 between somatic and normal samples. Tumor cellularity and purity estimates were determined using Sequenza and ABSOLUTE (25, 26).

HLA class I calling was performed on normal tissue whole exome sequencing from TCGA using OptiType, which has been shown to significantly outperform first-generation HLA prediction tools (27, 28) and mapped using RazerS3 as described previously (28, 29). HLA class II calling was performed using HLAreporter (30). Samples with low-certainty typing results were excluded. To identify potential novel peptides expressed in tumor tissue, DNA variants for each disease type were filtered for missense mutations in genes with a normalized expression count (scaled RSEM raw count) of greater than 1. After filtering, a sliding window method was used to identify 8–15-mer peptides centered on the mutation site. Peptides were then ranked for binding to each of the patients' specific HLA class I (8–14-mers) or II alleles (15-mers) using the consensus method provided by the Immune Epitope Database and Analysis Resource (31). Peptides with a median half-maximum inhibitory concentration of less than 50 nmol/L (class I) or a percentile rank of less than 1 (class II) were classified as potential neoepitopes.

Using publicly available data from 134 primary tumor resection samples, we profiled the genomic and transcriptional landscape of human PDA in the context of the immune microenvironment. We focused on a validated gene expression signature of granzyme A (GZMA) and perforin-1 (PRF1) to assess intratumoral cytolytic T-cell activity (cytolytic index; CYT; ref. 21). GZMA is a tryptase that induces caspase-independent programmed cell death and PRF1 is a pore-forming enzyme that mediates entry of granzymes into target cells, both produced by activated cytolytic CD8⁺ T cells and upregulated following response to immunotherapy (32–35). Although Rooney and colleagues pioneered the utility of this cytolytic index broadly across many human cancers in their initial report, pancreatic cancer was not available (21). Here, to assess cytolytic index in PDA, we obtained RNA sequencing data from The Cancer Genome Atlas (TCGA) for multiple tumor types including recently released data for PDA samples (36). The exact TCGA PDA samples included in our analysis are noted in Supplementary Table S1. Consistent with previous findings, we found that cytolytic index was highest in kidney, lung, and cervical cancers and lowest in glioblastoma, ovarian, and prostate adenocarcinoma (Fig. 1A; ref. 21). The median cytolytic index of PDA samples was comparable with that of other cancer types, including lung squamous cell carcinoma and stomach adenocarcinoma (Fig. 1A) and cytolytic activity in PDA was significantly higher than activity in normal pancreas (Fig. 1A, inset). Interestingly, while PDA has a median cytolytic index similar to stomach adenocarcinoma (8.59 ± 7.5 vs. 8.01 ± 11.9 CYT index), the distribution of cytolytic index is significantly narrower in PDA (Kolmogorov–Smirnov test, Supplementary Fig. S1).

To classify the subpopulations of PDA tumors with high versus low cytolytic activity, we stratified the PDA dataset by defining tumors in the top 10th percentile by cytolytic index as cytolytic-high (CYT High) and tumors in the bottom 25th percentile as cytolytic-low (CYT Low; Fig. 1B). Cytolytic-high PDA tumors were enriched for gene sets associated with activated CD8⁺, PD1^high T cells (37, 38), confirming that expression of GZMA and PRF1 correlated with immune response and infiltration of CD8⁺ cytolytic T cells (Fig. 1C; Supplementary Figs. S2 and S3).

As these data suggest that stratification based on cytolytic T-cell infiltration, as measured by the cytolytic index, may be associated with distinct PDA subtypes, we determined whether cytolytic activity is associated with genomic and transcriptional metrics of PDA tumor biology. Recent studies from multiple groups have demonstrated the extensive genetic and transcriptional diversity of PDA tumors (6–10). PDA can be stratified into at least three tumor subtypes based on gene expression profiling: (i) classical/pancreatic progenitor; (ii) squamous/quasi-mesenchymal/basal-like; and (iii) ADEX (aberrantly differentiated endocrine exocrine)/Exocrine-like (6–8). Furthermore, these PDA tumors types can overlap with gene programs associated with distinct stroma populations: (i) activated stroma, (ii) normal stroma, and (iii) immune gene signatures (6, 8). We assessed enrichment of gene programs defining PDA subtypes using the TCGA PDA dataset (Fig. 2A) and their association with cytolytic index (Fig. 2B). Cytolytic-high tumors were statistically enriched for the immune gene programs (GP7 and GP8) from Bailey and colleagues and the normal stroma gene program from Moffitt and colleagues (Fig. 2B; refs. 6, 8). Immune gene programs GP7 and GP8 contain markers for macrophages and T-cell coinhibition (GP7) and CD8⁺ T cells and B cells (GP8; ref. 6), and the normal stroma gene signature contains markers of pancreatic stellate cells (8). Pancreatic stellate cells have been linked to an immunosuppressive tumor microenvironment through PTX3 regulation of immune escape by blocking antigen presentation (39). Expression of genes defining the normal stroma gene program was increased in cytolytic-high tumors, suggesting a relationship between stromal microenvironment and T-cell infiltration (Fig. 2C). Cytolytic-low PDA tumors had statistical enrichment of gene programs associated with the classical/pancreatic progenitor tumor subtypes (Fig. 2B). Genes involved in pancreatic differentiation were increased in cytolytic-low tumors suggesting an inverse relationship between differentiation status and immune reactivity (Fig. 2D). Moreover, using a previously defined “immunome” gene signature of 28 distinct immune-cell–specific markers (40), we found that cytolytic-high tumors were associated with multiple other immune cell signatures (Supplementary Fig. S4). These data suggest that stratification of patients with PDA based on transcriptional profiling can differentiate between tumors with strong cytolytic T-cell response and tumors for which a privileged immune microenvironment precludes such responses (12).

We next sought to determine whether cytolytic activity correlated with distinct mutational profiles characterized for PDA (9, 10, 40). Curated mutational data for PDA was obtained from TCGA and we identified significantly mutated genes occurring in cytolytic-high and -low PDA tumors. Cytolytic-high tumors had a statistically significant association with mutations in TGFβR1/TGFβR2 as well as HMGB3 (Fig. 3A). Cytolytic-low tumors were associated with non-silent mutations in CDKN2A, ANKRD36, NCOA3, and HIST1H1B. Most mutations across the dataset were associated with G>A and C>T transitions, and the frequency of specific substitutions did not differ between cytolytic-high and -low tumors (Fig. 3B). Likewise, there was no association between cytolytic index and KRAS mutation type (Fig. 3C).

PDA is characterized by increased genomic instability with extensive copy number alterations in both human patients and genetically engineered mouse models (9, 10, 41). We next obtained the GISTIC2.0 analysis for PDA from TCGA and assessed copy number alterations between cytolytic subtypes. Cytolytic-low (but not high) tumors had recurrent copy number alterations at loci important in PDA (10), including MYC, NOTCH2, and FGFR1 (Fig. 3D, right). Consistent with increased copy number alteration and MYC amplification, cytolytic-low tumors had increased expression of gene signatures associated with genomic instability (42) and MYC target genes (43; Supplementary Fig. S5). Recurrent deletions were observed in cytolytic-low tumors at loci containing CDKN2A/B and SMAD4. MYC amplification has been observed in mouse models of hepatocellular carcinoma and human melanoma tumors and associated with reduced T-cell infiltration and cytolytic activity, suggesting that genomic events may modulate inflammatory response in PDA (44, 45). Cytolytic-high PDA tumors did not have recurrent copy number losses but rather had amplifications at 4q13.1 (TECRL), 9p13.3 (CA9, TPM2, C9orf100), and 18q11.2 (IMPACT, OSBPL1A; Fig. 3D, left). Assessed globally, cytolytic-high tumors had significantly fewer somatic copy number alterations (SCNA events) than cytolytic-low tumors (Fig. 3E and F).

Mutational analysis of tumor samples can be hampered due to tumor cellularity because this reduces the ability to confidently detect somatic mutation and copy number alterations. To ensure that the TCGA PAAD cohort has sufficient cellularity for our mutational analysis, we calculated cellularity estimates using the standard method ABSOLUTE (26) and confirmed using a second method Sequenza (25), although recognizing that these bioinformatics approaches are being continually refined. Tumor cellularity estimates did not correlate with either total mutation load or total copy number events in the TCGA PAAD cohort (Supplementary Fig. S6). Furthermore, TCGA PAAD cellularity estimates (59% ± 16.9%) were not significantly different than lung adenocarcinoma (LUAD, 52% ± 11.8%) and stomach adenocarcinoma (STAD, 53 ± 12.8%) tumor cohorts (Supplementary Fig. S7). Moreover, there was no difference in cellularity estimates between cytolytic subtypes, suggesting that observed differences in copy number and mutational load were not a likely result of variable tumor cellularity (Supplementary Fig. S8). Thus, distinct mutational and structural changes in the genome distinguish those PDAs with low versus high cytolytic activity.

We then determined whether cytolytic activity correlated with neoepitope load in PDA, as has been widely suggested for cancers in general (21). Neoepitopes, derived from peptides encoded by somatic tumor mutations and thus not subject to central tolerance in the thymus, have been demonstrated to preferentially drive T-cell recognition of tumor cells (46, 47). To determine whether cytolytic activity is associated with the presence of neoepitopes, we determined the frequency of total missense mutations and predicted those with potential to function as T-cell neoepitopes across the PDA TCGA dataset. To do this, we established a bioinformatics pipeline that predicts the ability of all possible 8–14-mer (HLA class I) or 15-mer (HLA class II) peptides inclusive of a tumor somatic missense mutation to bind to patient-specific MHC molecules with high affinity (<50 nmol/L for class I, <1% rank for class II). Peptides were then filtered by expression using TCGA RNAseq data (>1 normalized expression count). A total of 1.1 × 10⁴ unique variants leading to 5.1 × 10⁶ potential peptides were evaluated. This analysis revealed no correlation between total mutations per individual tumor and cytolytic index, although when viewed as a group, the cytolytic-low tumors exhibited a slight increase in the number of mutations per tumor compared with cytolytic-high tumors (Fig. 4A). Neoepitope load did not correlate with cytolytic activity in PDA, with striking examples of cytolytic-low tumors with multiple predicted neoepitopes and cytolytic-high tumors with few neoepitopes (Fig. 4B–E). Consistent with the findings from the overall mutation rate, the tumor cellularity estimates did not correlate with the number of MHC class I or II neoepitopes in PDA (Supplementary Fig. S6). In contrast, as expected (21), both lung adenocarcinoma and stomach adenocarcinoma showed a strong correlation between total number of mutations, predicted MHC class I neoepitope load, and number of mutations generating one or more neoepitopes in cytolytic-high tumors (top 10th percentile vs. bottom 25th percentile ranked by cytolytic activity; Supplementary Fig. S9). Stratification of patients with PDA based on established gene signatures (6, 8) also did not associate with increased neoepitope load in any PDA transcriptional subtype (Supplementary Fig. S10). Taken together, these data suggest that cytolytic activity in PDA, in contrast to other tumor types, is not driven by increased mutation or neoepitope load.

The tumor microenvironment in PDA contains a rich cytokine milieu with both pro- and anti-inflammatory factors that can regulate tumorigenesis (48). We therefore hypothesized that the expression of these cytokines and chemokines would be increased in cytolytic-high tumors across all TCGA tumor types. Consistent with this notion, we found that the expression of numerous pro- and anti-inflammatory cytokines and immune checkpoint molecules was significantly increased in cytolytic-high tumors across TCGA, including PDA (Fig. 4F). Specifically, cytokines previously shown to correlate with cytolytic index, including C1QA, C1QB, C1QC, CXCL10, and CXCL9, were differentially expressed in cytolytic-high TCGA tumors, including cytolytic-high PDA (21). The expression of regulatory T cell (Treg) markers was also significantly increased in cytolytic-high PDA tumors (Fig. 4I).

Finally, we assessed whether cytolytic-high PDA tumors exhibit increased expression of immune checkpoint pathways. We created an inhibitory checkpoint index to assess the expression of key checkpoint molecules across patients with PDA (Fig. 4G). The inhibitory checkpoint index was generated by taking the log-average expression in TPM of the following molecules: ADORA2A (A2AR), CD274 (PD-L1), PDCD1 (PD1), CTLA4, HAVCR2 (TIM3), IDO1, IDO2, PDCD1LG2 (PD-L2), TIGIT, VISTA (C10orf54), and VTCN1 (B7-H4). Using this index, we found a strong correlation between cytolytic activity and the expression of inhibitory checkpoint genes in patients with PDA, suggesting that as in melanoma (49), immune response in cytolytic-high tumors elicits multiple host and tumor mechanisms of immune suppression in the tumor microenvironment (Fig. 4H). However, CD274 (PD-L1) expression was uniformly low in PDA and was not differentially expressed between cytolytic subsets in PDA despite being increased in cytolytic-high tumors in the other TCGA datasets examined (Fig. 4F arrow). Multiple other inhibitory checkpoint molecules were expressed at markedly higher levels in cytolytic-high tumors (Fig. 4J). While PD-L1 expression was not changed between cytolytic subsets in PDA, the expression of other immune checkpoint molecules, including IDO1, IDO2, CTLA4, and PD-L2, were differentially expressed in cytolytic-high PDA tumors (Fig. 4J).

When compared individually, the expression levels of PDCD1, PDCD1LG2, CTLA4, IDO2, A2AR, TIGIT, and LAG3 have the highest correlation with the cytolytic index and, conversely, PD-L1 and VTCN1 (B7-H4) have the lowest correlation (Supplementary Fig. S11). These data suggest that multiple immune checkpoint pathways, other than the PD-1 axis, may mediate peripheral tolerance and immune escape in PDA; combinatorial targeting of these pathways may expand clinical benefit for patients with PDA.

Sophisticated approaches continue to provide unprecedented resolution of the immunobiology of human cancer. In this study, we performed an extensive integrated analysis of the transcriptional and genetic landscape of PDA in the context of cytolytic immune activity. By stratifying patients with PDA based on a validated cytolytic gene expression signature (not previously applied to PDA), we identified a small subset of patients with evidence of prominent T-cell reactivity. Beyond strong associations between cytolytic index and recently established transcriptional and genetic subtypes of this disease, our analysis revealed that low cytolytic activity tracked with increased genomic structural variation, most notably prominent and recurrent MYC amplifications and non-silent mutations and/or deletions in CDKN2A/B. Other distinct chromosomal aberrations were associated with cytolytic high PDA tumors. These data point to an underappreciated link between genomic alteration and immune activation in PDA, suggesting that genomic structural variation implicated in PDA progression may also fundamentally influence de novo or therapeutic antitumor immune activation, independently of host immune factors.

We also report the first characterization of neoepitope load in PDA, finding in a large sample size encompassing the entire TCGA dataset for PDA that high cytolytic activity failed to correlate with increased load of nonsynonymous mutations or predicted neoepitopes. This is in striking contrast to the correlation between cytolytic index and mutational burden in other tumor types such as lung and stomach adenocarcinoma (21), in which PD1 antibodies trigger clinically significant tumor regression. Indeed, recent studies have highlighted the primacy of tumor neoepitopes in T-cell recognition of tumor cells (46), renewing interest in patient-specific approaches such as personalized vaccines (50). Accordingly, the emerging paradigm is that neoepitope load determines sensitivity to immune checkpoint blockade (51–53). The lack of an association between neoepitope load and cytolytic index in PDA may reflect a tumor immunobiology that is distinct from that present in checkpoint blockade–sensitive tumors. Consequently, the assumption that neoepitope load is invariably associated with greater adaptive immunity may need to be reassessed, especially as it applies to patient selection in future PDA clinical trials. This lack of association in PDA may be due to extreme Kras-driven immunosuppression resulting in neoepitope-specific T cells that are not triggered, fail to expand, fail to infiltrate the tumor, or all of the above (54, 55). Alternatively, the inherently low mutation rate seen in PDA, compared, for example, with lung and stomach adenocarcinoma, may contribute to the lack of association between cytolytic reactivity and neoepitope load. These insights may explain the lack of clinical response to single-agent anti-PD1 therapy in patients with PDA and also indicate a need to vaccinate against tumor antigens and generate robust antitumor T cells to sensitize patients to checkpoint blockade therapy. However, if intrinsic tumor suppression can be overcome by these combinatorial methods, potent neoepitope-directed responses may yet be possible in the majority of PDA patients whose tumors contain dozens or more of these potential high affinity targets.

Our findings suggest that intrinsic oncogenic processes, rather than the availability of favorable immune targets, perhaps due to low mutation rate in PDA, may be the primary driver of immune activity in human PDA. MYC amplification is associated with decreased T-cell infiltration in mouse models of hepatocellular carcinoma and human melanoma tumors with low cytolytic activity, suggesting that genomic alterations besides neoepitopes can modulate inflammatory response (44, 45). Our comparison of the cytolytic index to established signatures of distinct PDA tumor and stroma types revealed a strong association between tumors with high cytolytic activity and a gene signature program representing a “normal stroma” phenotype (8), suggesting that the stromal microenvironment in PDA plays an important role in modulating inflammatory response. In PDA, oncogenic processes other than mutation and neoepitope load (e.g., recurrent MYC/NOTCH amplification) appear to be associated with a state of immune privilege that precludes host T-cell infiltration. Thus, immune checkpoint molecule upregulation in PDA may not be evidence of preceding T-cell immunity, as is likely the case in melanoma and lung adenocarcinoma.

Finally, we found that tumors with high cytolytic activity exhibited increased expression of multiple immune checkpoint genes such as CTLA4, TIGIT, TIM3, and VISTA. In contrast, PD-L1 expression was uniformly low. Whether therapeutic vaccines can generate sufficient cytokine-producing T cells that infiltrate PDA tumors and upregulate PD-L1 remains to be seen in clinical trials; regardless, our data provide a rationale for prioritizing immune checkpoints other than only PD-L1/PD-1 as therapeutic targets. While redundant suppressive mechanisms may undermine efforts to target any single immune checkpoint pathway, it may be possible to use expression profiling on a patient-specific basis to accomplish “immune precision medicine.” In addition, the lack of immune activation in response to high neoepitope load additionally suggests that immune checkpoint blockade may be insufficient in a relatively low-mutation rate tumor type such as PDA, and could be more effectively paired with immune-activating strategies that have shown promise is murine models (56, 57).

In summary, our findings suggest that it will be important to look beyond standard neoepitope-based strategies for immunotherapy in PDA and to focus instead on other tumor-intrinsic features that render these tumors immune privileged (12). The extent to which this immunobiology characteristic of human PDA also manifests in other types of carcinoma requires further investigation.

No potential conflicts of interest were disclosed.

The results shown here are in part based upon data generated by the TCGA Research Network (http://cancergenome.nih.gov/).

Conception and design: D. Balli, A.J. Rech, B.Z. Stanger, R.H. Vonderheide

Development of methodology: D. Balli, A.J. Rech, R.H. Vonderheide

Acquisition of data (provided animals, acquired and managed patients, provided facilities, etc.): D. Balli, A.J. Rech

Analysis and interpretation of data (e.g., statistical analysis, biostatistics, computational analysis): D. Balli, A.J. Rech, B.Z. Stanger, R.H. Vonderheide

Writing, review, and/or revision of the manuscript: D. Balli, A.J. Rech, B.Z. Stanger, R.H. Vonderheide

Administrative, technical, or material support (i.e., reporting or organizing data, constructing databases): D. Balli

Study supervision: B.Z. Stanger, R.H. Vonderheide

We are grateful to Anand Srinivasan, Jason Hughes, and members of the Stanger and Vonderheide laboratories for feedback and assistance.

This work was supported by the NIH (R01-CA169123; to B.Z.S. and R.H.V.), the American Cancer Society (postdoctoral fellowship to D.B.), and the Parker Institute for Cancer Immunotherapy (to A.J.R. and R.H.V.). Computational analysis was conducted on the Penn Medicine Academic Computing Services High-Performance Computing (Penn HPC) environment with the support of NIH grant 1S10OD012312.

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.