Genome-scale Pathway Flux Analysis Predicts Efficacy of anti-PD-1 Therapy in Melanoma

Background : Currently, predicting treatment efficacy to immunotherapy has been under extensive investigation. However, a putative biomarker for immunotherapy in melanoma remains to be found. Methods : Utilizing genetic data from two independent melanoma patient cohorts treated with anti-PD-1 therapy, the study described herein conducted a genome-scale pathway flux analysis (GPFA) and related statistical methods to determine whether specific pathways could be identified that are relevant to immunotherapy efficacy. Results : The analysis results highlighted three mechanisms responsible for the efficacy of immunotherapy in melanoma including 1) proper cellular functions in immune cells; 2) angiogenesis for the development and differentiation of immune cells; 3) energy metabolic remodeling to meet the activation of immune cells. Based on these discoveries, a pathway flux-based biomarker (IM.Index) was developed and assessed to validate its predictive ability with odds ratio (OR) of 3.14 (95%CI: 1.16-8.45; p=3.10E-3), sensitivity 76% and specificity 89%. The IM.Index achieved an objective response rate (ORR) of 70%. Comparison to other four putative biomarkers (TMB, NAL, neo-peptide load and cytolytic score) showed a comparative outcome with an hazard ratio (HR) of 1.83 (95%CI: 1.26-2.67; p=1.62E-3) and area under curve (AUC) of 0.82. Conclusion : These results indicate a translational potential of IM.Index, as a biomarker, for anti-PD-1 therapy in melanoma and the GPFA might pave a new path for biomarker discovery in immunotherapy. “collaborate” in an inseparable manner to ensure the functionality of the immune system 38 . A genome-scale pathway analysis might therefore provide insight into the functional status of the immune system and elucidate key mechanisms specific for anti-PD-1 therapy efficacy. This study was conducted in order to face these challenges utilizing genetic data from two independent melanoma patient cohorts treated with anti-PD-1 therapy (nivolumab & pembrolizumab) 33 34 . A genome-scale pathway flux analysis was developed and used to analyze the pathway activities on the data from the two cohorts, which investigated a new perspective in determining treatment efficacy on a genome-scale pathway level.

During cancer development diverse immune evasion mechanisms are employed to skip, prevent, or even terminate the host immunity so that cancer cells enter into the immune escape phase 9 . Based on the latest immunological research, anti-PD-1 therapy aims to inhibit the function of the B7-H1/PD-1 pathway implicated in diverse cellular functions including apoptosis 10 , angiogenesis 11 , proliferation 12 , anergy 12 and others 13 14 . Thus, anti-PD-1 therapy could impair all types of immune cells including Tcell 14 , B-cell 15 , dentritic cell 16 , marcophage 17 and NK cell 18 19 . A number of studies were conducted to determine the molecular and pathological mechanisms responsible for the efficacy of anti-PD-1 therapy. For instance, many studies focused on the leukocyte count within TME such as tumorinfiltrating leukocytes (TILs) 20 , tumor-associated neutrophiles (TANs) 21 , and tumor-associated macrophages 22 . Unfortunately, this type of quantification elicited contradicting results regarding the prognosis of immunotherapy 23 24 25 . More recent studies proposed that a high number of clonal neoantigens or mutational load may be associated to immunotherapy response. Tumors of this kind may harbor a great number of mutated immunogenic neoantigens that are able to be recognized through the host immunity 26 27 . This finding led to diverse biomarker-related studies in mutational burden (TMB) 28 and neoantigen load (NAL) 29 , or TCR repertoire profiles 30 and assays of checkpoint status 31 32 . However, a clear indication or a putative functional biomarker for anti-PD-1 therapy remains to be found. Interestingly, two recent studies generated genome-scale data of melanoma patients treated with anti-PD-1therapies and analyzed the possible associations of treatment efficacy on a DNA and RNA level 33 34 . However, possible mechanisms for the efficacy of immunotherapy in melanoma could not be singled-out.
What most of the above mentioned studies have however failed to consider is that any activation of immune system is an energy demanding process interconnected through diverse intensive signaling procedures 35 36 37 . Thus, both metabolism and signaling "collaborate" in an inseparable manner to ensure the functionality of the immune system 38 . A genome-scale pathway analysis might therefore provide insight into the functional status of the immune system and elucidate key mechanisms specific for anti-PD-1 therapy efficacy. This study was conducted in order to face these challenges utilizing genetic data from two independent melanoma patient cohorts treated with anti-PD-1 therapy (nivolumab & pembrolizumab) 33 34 . A genome-scale pathway flux analysis was developed and used to analyze the pathway activities on the data from the two cohorts, which investigated a new perspective in determining treatment efficacy on a genome-scale pathway level.

Study Design and Objective
This retrospective study focused on the genome-scale pathway analysis of melanoma patients treated with anti-PD-1 therapy. We selected nivolumab and pembrolizumab as a representation of anti-PD-1 therapy. Genetic data of two previously published study cohorts of melanoma patients were applied for the analysis purpose. The objective of this study was to identify common molecular mechanism of nonresponder in anti-PD-1 therapy and summarize these mechanisms into a validated biomarker. Figure 1C visualizes the workflow of the intended study. The need for ethics approval was waived because the published data was applied.

RNAseq of Patients and Data Preprocess
This study utilized two cohorts of melanoma patients from two previously published independent studies: Riaz et al., 2017 (the cohort C1) 33 and Hugo et al., 2016 (the cohort C2) 34 . Both sets of RNAseq were downloaded from the GEO (www.ncbi.nlm.nih.gov/geo/) with corresponding access IDs: GSE91061 and GSE78220. For the preprocess of the RNAseq, sequence reads were aligned using Salmon (v0.9.1), and the transcript quantification was performed by aggregating the transcript counts.
Salmon uses the raw fastq sequence reads and the reference transcriptiome assembly (GRCh37). A quasi-mapping method is then applied where the programm tool computes the mapping of reads to the transcript positions without performing a base-to-base alignment of reads to the transcript.

Genome-Scale Pathway Flux Analysis (GPFA)
The GPFA is developed to detect whether a therapeutic intervention can invoke a genome-scale significant change of data flux within a cellular system. The data flow (flux) was generated from AutoAnalyze 39 based on related readout components in an incorporated molecular model, a modified version of the molecular signaling map (MSM) 40 . The model MSM possesses four layers with following relationship: (1) gene→ (2) RNA→ (3) protein/complex/compound → (4) pathway → (1) gene (Supplement Figure 2). The last layer (78 pathways) is connected back to the first layer (2938 gene) due to diverse feedback mechanisms including transcriptional and translational regulations. For instance, a molecular reaction R with role ∈ {e, g, i, s, tr(a) , tr(r) }, R symbolizes a reaction set; (e: enzyme; g: gene; i: inhibitor; s: substrate; tr(a): transcriptional activator , tr(r): transcriptional repressor): Ij (reaction,role) = ∏ objects ∈ reactant (reaction,role).objects (1) c(object) The data-flow for this reaction I ∈ R is computed by applying the mass action law with the required input concentrations of related reactants, c(objects). Here, c(objecti) = ∏ reactants ∈ (reaction,role).genes (2) c(gene) The concentration of each reactant is determined by the input of gene expression profile of a corresponding patient. The flux of a pathway P is defined as flux(P), the amount of information flowing through this pathway: flux(P) = ( ∑Ii(reaction, role) / N(P) ) -flux( Crosstalk(P) ) i ∈ P N(P): the number of reactions of the pathway P; Crosstalk(P): the crosstalks pathways related to the pathway P.

The Definition of IM.Index
The IM.Index was defined as a score index to summarize the flux from multiple immune system related pathways and its application was to predict the resistance of anti-PD-1 therapy. A network-oriented data-flow analysis, AutoAnalysis 39 , was applied to conduct the calculation of IM.Index. The formula of an IM-Index of a given patient with gene expression profile (GEP) is defined following (Supplement Information): p ∈ signaling transduction p ∈ energy metabolism

Optimization of Parameters of IM.Index
The i-th GEP from the responder group and j-th from the non-responder group are GEPi and GEPj respectively. The difference of IM.Index (Diff) between both groups is defined as following: i ∈ responder group j ∈ non-responder group Diff(i, j) reaches its maximum value for the estimation of parameters, α and β. We applied stochastic gradient descent method (SGD) to solve the optimization issue of (5) and used five-fold crossvalidation to avoid the parameter overfitting. For the optimization, α and β were set to 1.78E-4 and 2.37E-4.

Statistical analysis
Benferoni-correction was applied to define whether the changes of pathway flux are statistically significant between different groups including CR, PR, SD and PD from both cohorts. Wilcoxon test was applied to calculate the statistical significance whether IM.Index could differentiate between the responder and non-responder group. The prognostic impact of classifier was evaluated by the Kaplan-Meier method and the log-rank test. Receiver operating characteristic (ROC) analyses were then done to compare their performance. Area under the curve (AUC) and the 95% confidence intervals were calculated. With the null hypothesis that the biomarker gives equal or better classification than the index, one-sided hypothesis test was performed and bootstrapped p-value of the test was calculated (with 2,000 replications). Studenten t-test was performed to investigate whether the changes of IM.Index in melanoma between before and after therapy was significant. Multivariable cox model were applied to assess the association of different potential biomarkers with overall survival (OS). For pathway flux analyses and all other analysis, p values <0.1 and p values <0.05 were considered to be statistically significant respectively.

Patients
The first cohort of melanoma patients (n=42) were treated with nivolumab in accordance with the protocol 33 . Treatment results were divided into four categories: complete response (CR, n=3), partial response (PR, n=6), stable diseases (SD, n=15), and progressive disease (PD, n=18). The gene expression data of this cohort was collected before and after treatment. This cohort is referred to as "C1" during this study. The second set of melanoma patients (n=28) were treated with pembrolizumab in accordance with the treatment protocol 34 . Treatment results were divided into three categories: complete response (CR, n=5), partial response (PR, n=10) and progressive disease (PD, n=13). The gene expression data of this cohort was collected before the treatment. This cohort is henceforth referred to as "C2". For this retrospective study, CR and PR were defined as a responder group, whereas SD and PD (or PD alone) were defined as a non-responder group. The patient demographics can be found in the original study publications. Both were similar in regard to sex, age and cancer stage distribution.

Molecular Pathway Difference Between Responder and Non-Responder in Nivolumab-treating Melanoma
In order to determine whether molecular pathway mechanisms exist between responder and nonresponder of anti-PD-1 therapy in cohort C1, the GPFA was conducted between the responder and nonresponder group in C1 before therapy (pre-therapy). The result showed that multiple signaling pathways including death-receptor, insulin-receptor, NFAT, RET, TCR, ErbB, VEGF, and JAK-STAT displayed a significantly higher flux in the responder group before therapy (p<0.1; Fig. 1A; Table 1).
Among them, signaling death-receptor, NFAT and VEGF were associated with cellular processes such as apoptosis 41 , cell cycle 42 and angiogeneiss 43 . Insulin-receptor, ErbB, VEGF, JAK-STAT and TCR on the other hand were highly related to immune cell functionality including T-cell 44 45 and B-cell 46 47 . Furthermore, one energy-relevant metabolic pathway and fatty.acids showed a significantly higher flux in the non-responder group before therapy, whereas pantothenate.CoA displayed a higher flux in the responder group (p<0.1; Fig. 1A; Table 1). After therapy (on-therapy) GPFA was compared between the responder and non-responder group in cohort C1. The results showed that these aforementioned signaling pathways continued to show a high activity status in the responder group and additionally in several patients in SD group, with the exception of JAK-STAT and death-receptor pathways (Fig. 1B). Pantothenate.CoA also remained an increased flux in the responder group (p<0.1; Fig. 1B; Supplement Info.). Interestingly, two metabolic pathways ether.lipid and glycerolipid displayed a high flux in the responder group after therapy (p<0.1; Fig. 1B; Table 1)

Molecular Pathway Difference Between Pre-and On-Therapy in Nivolumab-treating Melanoma
A refined analysis focused on the functional alternations of molecular pathways in cohort C1 before (pre) and after (on) anti-PD-1 therapy. At first, GPFA analysis was performed in following four subgroups separately: CR_pre vs. CR_on ( Fig. 2A), PR_pre vs. PR_on, (Fig. 2B) SD_pre vs. SD_on (Fig.  2C), and PD_pre vs. PD_on (Fig. 2D). The results showed that in the CR subgroup two signaling pathways and four metabolic pathways showed a significantly increased flux following anti-PD-1 therapy (p<0.1; Fig. 2A; Table 2). It is noteworthy that all these six pathways are essential to the function of T-cells 48 49 50 . The comparison between pre-and on-therapy in PR subgroup revealed that multiple signaling and metabolic pathways were upregulated after therapy (p<0.1; Fig. 2B; Table 2). The six pathways reported in Fig. 2A also appeared in Fig. 2B. The result in non-responder group (SD+PD) showed that only diverse metabolic pathways displayed a higher flux after therapy and none of signaling pathways showed a significanly changed flux ( Fig. 2C & 2D). For instance, fluxes of several Vitamin B-related metabolic pathways (Folate, Riboflavin, Thiamine) were significantly altered after therapy (p<0.05; Fig. 2A; Table 2). Among them, thiamine metabolism was only upregulated in this group, not in responder group after therapy. Interestingly, metabolic resource-relevant pathways nicotinate, glutamine, arginine, and mannose showed a high flux only in the non-responder group after therapy, too (p<0.1; Fig. 2A; Table 2). Furthermore, fluxes of pentose-phosphate pathway and fatty acid metabolism were significantly reduced only in the non-responder group after therapy (p<0.1; Fig. 2A; Table 2). Fig. 2E summarizes the results of these four sub-figures (2A-2D).

Molecular Pathway Difference Between Responder and Non-Responder in Pembrolizumab in Melanoma
Here we investigated the functional alternations of molecular pathways in the cohort C2. Before pembrolizumab treatment GPFA was conducted to compare responder and non-responder in the cohort C2. The results showed that multiple signaling pathways had a higher flux activity in non-responder group (Fig. 3A). For instance, the immune system-based pathways IFN, TLR and JAK.STAT (p<0.05; Table 2); proliferation-related pathways mTor, MAPK, WNT (p<0.05; Table 2); apoptosis-related pathways TGF-beta and p53 (p<0.05; Table 2); angiogenesis-related pathways PDGF, RAR and renin.angiotensin (p<0.05; Table 2). Interestingly, many metabolic pathways showed a increased flux activity in the non-responder group (Fig. 3B). For instance, diverse cellular nutrient-related metabolic pathways including purine, glutathione, tryptophan, arginine, thiamine (p<0.05; Table 2); several metabolic energetic pathways including citrate cycle, glycolysis, nicotinate, mannose and amino.sugar.and.nucleotide.sugar (p<0.05, Table 2).

Biomarker for anti-PD-1 Treatment of Melanoma Patients
Based on the pathway results of aforementioned GPFA an anti-PD-1 treatment-resistance score, termed IM.Index, was defined as a sum of flux from these immunotherapy-related signaling and metabolic pathways (Methods). Parameters of IM.Index were determined by maximizing the difference of IM.Index between the responder group and the non-responder group (Methods). IM.Index was calculated for each patient in both cohorts. Via IM.Index the responder-and non-responder groups in the cohort C2 were able to be differentiated with an odds ratio (OR) of 2.12 (95%CI: 1.22-3.66; p=4.50E-4; AUC=0.87; Fig. 4A-4B). Both groups with median of IM.Indexes were visualized in the pathway interaction map of leukocytes (PIML) respectively (Fig. 4C & 4D; Supplement Figure 1). The IM.Index predicted OS in this cohort with an hazard ratio (HR) of 1.24 (95%CI: 1.04-1.47; p=1.40E-2; Fig 4E). Subsequently the IM.Index was used to predict the response of anti-PD-1 therapy and OS in the cohort C1. The results showed that the IM.Index differentiated between the responder (median: 4.39) and nonresponder group (median: 5.41) in this cohort with an OR of 3.14 (95%CI:  Fig. 4L-4K). IM.Indexes (median) of both groups after therapy are visualized in the PIML (Fig. 4M & 4N). The prediction rate of OS by IM.Index remained significant with an HR of 1.35 (95%CI: 1.13-1.62; p=4.70E-2; Fig. 4O).

Comparsion Between IM.Index and Other Potential Biomarkers in PD-L1 Therapy
The comparison analysis between IM.Index and other four potential biomarkers (TMB, NAL, neo.peptide load, and cytolytic score) was performed according to terms of HR, ORR, and AUC for the cohort C1. Initially the cox's proportional hazard regression with multivariable model was performed investigate these five variables in the relation to OS of patients in the cohort C1. The results showed that IM.Index with an HR of 1.83 (95%CI: 1.26-2.67; p=1.62E-3; Concordance-Index=0.69) outperformed the other four biomarkers (Table 3; Fig. 5A). Moreover, TMB was showed to be related to OS (HR: 1.02; p=5.93E-2) and could be combined with IM.Index to improve the OS prediction in this cohort (Table 3). Subsequently the ORR analysis was performed for each of these variables. The result showed that IM.Index achieved an ORR of 70%, which is much higher than the currently observed ORR 40% to WHO criteria (Fig. 5B). As compared with other four biomarkers, the IM.Index outperformed them regarding ORR ( Fig. 5B; Table 3). Afterwards, response prediction for anti-PD-1 therapy was analyzed for all these variables. The results showed that IM.Index (AUC=0.82, 95%CI: 0.68-0.95) outperformed other four potential biomarkers regarding anti-PD-1 therapy in this cohort (Table 3 ; Fig. 5C). Interestingly, the cytolytic score (AUC=0.71) was shown as second to be predictive in the response of this therapy.

Discussions
Our study performed a genome-scale pathway flux analysis (GPFA) on the efficacy (Responder Vs. Non-responder) of anti-PD-1 therapy (Nivolumab & Pembrolizumab) in melanoma. Based on these results we were able to single-out three molecular mechanisms that could be responsible for the efficacy of anti-PD-1 therapy. Applying these mechanisms, we defined an anti-PD-1 treatment-resistance score (IM.Index) with which one could stratify individual melanoma patients who would or would not benefit from therapy (sensitivity 76%, specificity 89%). Furthermore, the IM.Index was significantly associated with OS in both cohorts of melanoma patients (p<0.05). The validation result showed that IM.Index outperformed other four potential biomarkers (TMB, NAL, neo.peptide load, and cytolytic score) for anti-PD-1 therapy in melanoma (Fig. 5).
Via a GPFA using the cohort C1 we identified several signaling pathways that showed a higher flux in the responder group, as compared to the non-responder group before and after anti-PD-1 therapy (nivolumab) (Fig. 1). Among these upregulated signaling pathways, death-receptor, NFAT, and VEGF are associated with cellular processes such as apoptosis 26 , cell cycle 27 and angiogeneiss 28 , whereas insulin-receptor, ErbB, VEGF, JAK-STAT, and TCR are highly related to immune cell functionality including T-cell 29 30 and B-cell 31 32 (p<0.1; Table 1; Fig. 1). This result may reflect the therapeutic effect of anti-PD-1 therapy (nivolumab), which clearly differentiates them from patients in nonresponder group. Furthermore, metabolic pathway pantothenate.CoA also showed a higher flux in the responder group before and after therapy (Fig. 1), which may suggest that coenzyme A (CoA) as a final product of this pathway, could serve as an essential indicator for the efficacy of anti-PD-1 therapy. In addition, we compared each of these sub-groups (CR, PR, SD, and PD) using the cohort C1 between pre-and on-therapy (Fig. 2). The results showed that in SD and PD subgroups only metabolic pathways such as glycolysis, valine.leucin.isoleucin-, arginine-, glutamine-metabolism and others were upregulated during on-therapy and no signaling pathways was affected. This may suggest that anti-PD-1 therapy (nivolumab) could result in irAEs in the non-responder group in this cohort instead of a therapeutic benefit. For the cohort C2 only genetic data before anti-PD-1 therapy (pembrolizumab) was collected 34 . The GPFA results of this cohort were divided into two parts: signaling-and metabolic perspective (Fig. 3A  & 3B). From a signaling perspective we observed that several pathways such as mTor, Death-receptor, IFN, TLR, MAPK, PDGF and other displayed a significantly higher flux in the non-responder compared to the responder group in this cohort (p<0.05; Table 2; Fig. 3A). Since several of these signaling pathways are intricately involved with the function and development of immune cells this result might propose a selection criteria for patients that would not respond to anti-PD-1 therapy (pembrolizumab) 36 51 52 . Some other pathways such as NFkB, RET, TGF-beta, WNT, NOTCH showed a moderately higher flux in the non-responder group (p<0.1; Table 2; Fig. 3A). It is noteworthy that several of them determine the function-duration and -strength of immune response 53 54 55 . Flux data from other pathways such as NFAT, TCR, and interleukin were similar between both groups (Fig. 3A). From a metabolic perspective we observed that multiple amino acid-based metabolic pathways such as alanine, cysteine, tryptophan, arginine, and others showed an increased flux in the non-responder group in this cohort, as compared with the responder group (p<0.1; Table 2; Fig. 3B). Furthermore, the amino acid, nucleotide and sugar interconversion pathway (amino.sugar.and.nucleotide.sugar, p<0.01; Fig.  3B) showed a consistently higher flux in this group. This result may strongly indicate that tumor cells in the non-responder group have already acquired a large amount of metabolic energy. In addition, N.glycan biosynthesis displayed a significantly higher flux in the non-responder group, this pathway is implicated in key pathophysiological mechanisms during T-cell development 50 56 . The thiamine pathway was also upregulated in the non-responder group, which may inhibit intracellular p53 activity to restrain its regulation on apoptosis 57 . Any deregulation of pathway activities may greatly contribute to the dysfunction of the immune system in the non-responder group, which are clearly distinct from those in the responder group to anti-PD-1 therapy (pembrolizumab).
Given the above analysis results, we summarized common pathway properties in both cohorts regarding nivolumab and pembrolizumab treatment and summarized potential molecular mechanisms that may be essential for a response of anti-PD-1 therapy. One of these pathway properties were the fluxes data from the pathways NFAT (p>0.1) and TCR (p>0.1), these were not significantly elevated in both non-responder groups ( Fig. 2E; Fig. 3A). This result may indicate that cellular processes in immune cells, including T-cells and B-cells, such as cell cycle 42 58 and apoptosis 42 59 could be interfered with. The VEGF pathway (p>0.1) did not display a high flux in both non-responder groups ( Fig. 2E;  Fig. 3A). This may indicate that insufficient angiogenesis conld impair the development and differentiation of T-cell 60 61 , and therefore impede the efficacy of both drugs in anti-PD-1 therapy. Flux of Insulin-receptor pathway (p>0.1) was not significantly enhanced in both non-responder groups ( Fig.  2E; Fig. 3A). This could imply that energy-demanding process fueled by glucose consumption during T-cell activation may be inhibited 44 , which may result in a significant reduction of the treatment efficacy.
Other observed pathway properties were that several energy-relevant metabolic pathways including glutamine, nicotinate, arginine, and mannose displayed a elevated flux in both non-responder groups (p<0.01; Fig. 2E & 3B). These may imply that tumor cells could still collect sufficient resources for proliferation 62 63 64 . Thiamine metabolism was upregulated in both non-responder groups (p<0.01; Fig.  2E & 3B). The upregulation of this pathway may directly inhibit the activities of p53 function 65 and other cellular relevant functions 66 , leading to the reduction of apoptosis in tumor cells. Pantothenate metabolism was only affected in the responder group in both treatments (p<0.01; Fig. 2E, 3B). This result may give an interesting hint into the responsiveness of anti-PD-1 therapy in the melanoma because this metabolic pathway is essential in the control of cellular CoA required for the proper function of the immune system 67 68 . Thus, the therapeutic effect of anti-PD-1 therapy (nivolumab & pembrolizumab) may require a sufficient availability of cellular CoA. In summary three molecular potential mechanisms for the successful application of anti-PD-1 therapy were identified through our analysis: 1) proper cellular functions in immune cells; 2) angiogenesis for the development and differentiation of immune cells; 3) energy metabolic remodeling and the availability of energy resources to meet the activation of immune cells.
Two previous studies (C1 and C2) separately analyzed the efficacy of nivolumab (C1) and pembrolizumab (C2) on a DNA and RNA level. However, our study performed pathway flux analysis on a genom-scale (GPFA). GPFA considers the relationships of four levels in biology: (1) gene→ (2) RNA→ (3) protein/complex/compound → (4) pathway → (1) gene. Diverse regulation mechanisms were considered within the GPFA, which includes transcriptional regulation and feedback control mechanisms controlling biological connections between (1) and (2) and between (4) and (1); translational regulation mechanisms controlling biological connections between (2) and (3). Thus, the application of GPFA offers a new perspective in drug efficiency research and its underlying therapeutic mechanisms on a genome-scale. For instance, GPFA revealed that none of signaling pathways were upregulated in the non-responder group in cohort C1 after therapy, which strongly and clearly indicates that nivolumab treatment could not provide any clinical benefits for this subgroup of patients (Fig. 2E).
Many recent studies explain that IFN pathway plays a key role in the resistance to immune checkpoint blockage therapy 69 70 71 . In line with these results, our study also showed that not only did this pathway display a high flux in the non-responder group (Fig. 3A), but also multiple other pathways such as VEGF, Insulin-receptor pathway and metabolic pathways (Fig. 3). This demonstrates the advantage of an investigation from a genome scale perspective. Isolated investigation of a single pathway may miss essential functional implications in the nature of cellular signaling transduction and crosstalk. Furthermore, diverse upregulated metabolic pathways in non-responder group from both treatments may also imply a high rate of irAEs to anti-PD-1 therapies. Therefore, the results from our study could fill the gap of these two previous studies 33 34 and the proposed IM.Index may reflect the strength of the host immunity. One major limitation of this study is the limited number of melanoma patients. Future studies should recruit more patients and conducted prospective studies to further verify the prediction ability of this approach.

Conclusion
Aberrant cellular signaling transduction and remodelling energy metabolism are fast becoming central factors in contributing to immunotherapy resistance. The pathway flux analysis on a genome-scale described above could provide a new prospective on the efficacy of anti-PD-1 therapy in melanoma patients. To our knowledge, our study could be the first one to introduce a system level pathway activity-based biomarker (IM.Index) for immunotherapy. However, the IM.Index should be further validated with multi-center data to improve and consolidate its clinical utility. We envisage that a similar approach could be developed to investigate the efficacy of immunotherapy for other types of cancer.

Ethics approval and consent to participate
The need for ethics approval was waived because the published data was applied in this study.

Consent for publication Not applicable
Availability of data and materials GEO (www.ncbi.nlm.nih.gov/geo/) with corresponding access IDs: GSE91061 and GSE78220

Competing interests
We clarify that there is no conflict of interest among co-authors to report.

Acknowledgment
We thank Dr. Yisong Wan and Dr. Michael Lauseker for their insightful discussion.       Overall survival analysis in the cohort C1 of melanoma patients between these 5 potential biomarkers. B: Comparison of objective response rate (ORR) between these 5 potential biomarkers. C: Comparison of AUC between these 5 potential biomarkers. Table 1. List of significant pathways during anti-PD-1 therapy (pre-and on-therapy) in melanoma. Table 2. List of significant pathways in comparison in complete remission, partial remission, stable disease and progressive disease during anti-PD-1 therapy (CR_pre VS. CR_on, PR_pre VS. PR_on, SD_pre VS. SD_on, PD_pre VS. PD_on, correspondingly) in the cohort C1 and between respondergroup and non-responder-group in the cohort C2. Table 3. Comparison of five potential biomarker for anti-PD-1 therapy in melanoma.