Hopf bifurcation in a CTL-inclusive HIV-1 infection model with two time delays

1.   Introduction

Head and neck squamous cell carcinoma (HNSCC) is a common cancer and develops from the oral cavity, oropharynx, larynx and hypopharynx ^[1]. There were 264,211 male and 113,502 female new cases arise from lip or oral cancer around the world during 2020 because there are no efficient and identified prediction and therapeutic approaches ^[2]. Traditional therapies are surgical resection, radiotherapy and chemotherapy, with many current anti-cancer and cancer metastasis agents achieving their therapeutic effect by promoting apoptosis of cancer cells. Therefore, exploring different forms of cancer cell apoptosis to resist cancer cell metastasis and finding more precise therapeutic sites may help improve the prognosis of HNSCC patients.

In 2015, researchers identified an apoptotic pathway induced by copper that alters intracellular ROS levels and mitochondrial inner membrane sites ^[3]. Subsequently, it was demonstrated that copper may induce apoptosis through specific pathways ^[4,5]. Peter Tsvetkov's team finally defined copper-induced cell death as cuproptosis, which resulted in lipoylated protein aggregation and subsequent iron-sulfur cluster protein loss in 2022 ^[6]. Copper could chelate diethyldithiocarbamate and induce oxidative stress-mediated cancer cell death, with copper ionophores preferentially inducing cancer cell proptosis due to their selective nature compared to normal cells ^[7]. The discovery of cuproptosis offers feasible therapeutic solutions for the targeting of cancers and anti-cancer agent resistance.

Long non-coding RNAs (lncRNA) acted as an essential regulatory function in gene transcription and recently have been shown to regulate cell death patterns such as apoptosis, autophagy, and ferroptosis and thus affect the prognosis of HNSCC patients. lncRNA ANRIL, lncRNA NEF lncRNA AC026166.2-001 can promote squamous cancer cell death in apoptosis, however, more lncRNA are upregulated and inhibit apoptosis in laryngeal squamous cell carcinoma tissues, such as lncRNA FTH1P3, lncRNA PVT1, lncRNA MALAT1 ^{[8,9,10,11,12,13]}. The down-regulation of lncRNA GACAT1 inhibited the proliferation and migration of oral squamous cell carcinoma and promoted cell apoptosis and autophagy ^[14]. LINC01207 silencing promoted apoptosis and autophagy in oral squamous cell carcinoma, and LINC01207 overexpression had the opposite result ^[15]. LncRNA GAS5 could activate autophagy in laryngeal squamous cell carcinoma, with enhanced autophagy-related proteins after GAS5 overexpression ^[16]. In addition, ferroptosis and cuproptosis are newly discovered patterns of metal ion-induced cell death ^[17]. Whereas models predicting the prognosis of HNSCC patients based on ferroptosis-related lncRNA have been constructed, the field of cuproptosis still exists indeed.

In this study, we explored the genomic profile of cuproptosis-related in cancer microenvironment through single-cell RNA sequence data from the GEO database. Based on the differentially expressed genes between cancer and normal patients from TCGA, we screened for cuproptosis-related lncRNA with the Pearson test. A cuproptosis-related lncRNA prognostic model was constructed in the training set and provided a model formula to score the risk for patients who were divided into high-risk and low-risk score groups. Enrichment functions, pathways, immunogenomic landscape analysis, and chemotherapeutic drug sensitivity analysis were explored for patients in high- and low-risk groups. Furthermore, the model combined clinical-pathological characteristics to screen the independent prognostic factors and constructed a nomogram. The workflow of this study was shown in Figure 1.

Figure 1.  Workflow of this study.

DownLoad: Full-Size Img PowerPoint

2.   Materials and methods

2.1.   Data extraction

We downloaded bulk RNA sequence data from the cancer genome atlas dataset (TCGA, https://portal.gdc.cancer.gov/) for 44 normal humans and 501 cancer patients and obtained overall survival and clinical-pathological characteristics for patients. The cancer single-cell RNA sequence data from GSE195832 were taken from the Gene Expression Omnibus (GEO, https://www.ncbi.nlm.nih.gov/geo/) and identified cuproptosis gene were extracted from previous research ^[6,18,19] (Table S1). The genes were confirmed to comprise the lipoic acid pathway, pyruvate dehydrogenase complex, protein lipoylation enzymes and copper transporter protein. In order to pinpoint more clinically valuable cuproptosis-related genes, we performed a differential analysis of the expression matrix of TCGA patients in the cancer and normal groups and intersected the results with cuproptosis-related genes. The 12 expressed cuproptosis-related genes were obtained and used as subsequent cuproptosis-related lncRNA screens (CDKN2A, GLS, NLPR3, GCSH, DLD, FDX1, LIAS, MTF1, ATP7A, DBT, DLAT, PDHB). We obtained the annotation file from the Genecode website (https://www.gencodegenes.org/human/) to distinguish lncRNA and mRNA in counts. The relationship between differentially expressed genes and lncRNA were analyzed through Pearson test using limma package. The 486 cuproptosis-related lncRNAs were extracted from the expression profile based on a coefficient value > 0.4 and P-value < 0.001. All sequencing data obtained from publicly available repositories were approved and monitored by local ethics committees.

2.2.   Single-cell RNA sequencing data processing

The four-cancer single-cell RNA sequencing data imported the Seurat package in R4.1.1 to exclude the low-quality cells of which total counts of UMI were less than 500 or more than 15,000 and more than 10% mitochondrial genes ^[20]. The normalization of highly variable genes was based on a negative binomial regression model with regularised parameters using the SCTransform function. The FindIntegrationAnchors and IntegrateData function were utilised to merge four normalisation sequencing datasets. Batch correction of the integrated data was performed using the RunPCA function in Seurat, and the datasets were projected into the 50 principal components. We calculated the spatial distances of the different cell subpopulations and projected them into two dimensions using the RunUMAP function. Each cell cluster in the UMAP dimensions was classified as a particular cell type (e.g., Cancer, CD8 T cell, CD4 T cell, Fibroblast) based on marker genes of the identified cell clusters from previous research ^[18].

2.3.   Enrichment function and pathway analysis

We performed the differently analysis between patients in the tumor and normal group, and between high-risk score patients and low-risk score patients. Differentially expressed genes (DEGs) of cuproptosis were screened by limma packages with an FDR < 0.05 to identify the valuable cuproptosis genes in HNSCC by limma package in R4.1 ^[21]. With ClusterProfiler's enrichGO function, DEGs were analysed by Gene Ontology (GO) enrichment function for three categories (cellular components, molecular functions, and biological processes) ^[22]. The enrichKEGG function in the ClusterProfiler package examined biological pathway of DEGs. The significant results of enrichment function and pathway analysis were selected with cutoff value of P < 0.05.

2.4.   Construction of the cuproptosis-related lncRNAs prognostic model

An analysis of 501 patients with OS and cuproptosis-related gene expression was performed on a separate set of 250 patients for training and 201 patients for testing. Univariate COX regression screened cuproptosis-related lncRNAs to determine the potential signature associated with prognosis with P-value < 0.05 in training set. Then, Lasso regression were applied to screen the cuproptosis-related prognostic lncRNAs signature using 1000 cross-validations, allowing for the screening of highly correlated genes and minimizing the risk of overfitting the prognostic model. Then, the regression coefficients of each cuproptosis-related prognostic lncRNA had been calculated. Finally, a cuproptosis-related lncRNAs prognostic model was constructed with 11 lncRNAs under optimal parameters through multivariate COX regression with P < 0.05. Based on the regression coefficients and the expression matrix, it could construct the risk score formula: Risk score = ${\sum }_{1}^{n}coefficients*Exp$, Exp means expression of each cuproptosis-related prognostic lncRNAs in FPKM format. The risk score was evaluated for HNSCC patients, and they were then assigned to high or low-risk groups based on whether it is higher than median risk score in each set. The OS and lncRNA expression data were imported into the survival and survminer packages and Kaplan-Meier curves were plotted in R4.1.1. With the support of the timeROC package, ROC curves of the model were plotted.

2.5.   Testing of the cuproptosis-related lncRNA prognostic model

The cuproptosis-related lncRNA prognostic model was tested simultaneously on the training, test and the entire sets. Scatter map were plotted to illustrate the distribution of risk score and survival status. Principal component analysis (PCA) and t Stochastic Neighbour Embedding (TSNE) analyses were conducted to show differentiation in two risk score patients with the Rtsne package in R software. The expression matrix of 11 lncRNA for patients of different risk groups were extracted for being normalized and scaled. We use the prcomp function to perform PCA analysis on the expression matrix according to different groupings. We putted the expression matrix into Rtsne function to perform TSNE analysis with perplexity as 10 and max_iter as 500. Kaplan-Meier analysis and ROC were performed for 1-, 3-, and 5-years in the all sets to examine the risk score. The package was utilized to process data, including survival, survminer, rms and timeROC packages.

2.6.   Construction of a predictive nomogram

The clinical-pathological characteristics and the model were are to screen independent prognostic signatures using univariate and multivariate COX regression. Moreover, the area under ROC curves was utilized to estimate diagnostic and prediction capability of model or clinical-pathological characteristics. The C-index (concordance index) evaluates the model by performing a pairwise analysis of clinical cases to assess the consistency of predicted OS with true prognostic outcomes. An RMS package is utilized to build a nomogram that predicts OS of HNSCC patients based on clinical-pathological traits and risk scores. The R package was utilised to process data such as rms, replot and survival packages.

2.7.   Immunogenomic landscape analyses and chemotherapeutic agent analysis

The immune functions and checkpoints between two risk score patient HNSCC patients were evaluated using ssGSEA algorithms to compute the bulk gene expression profile based on the GSVA package ^[23]. The cancer immune dysfunction and exclusion (TIDE) used to model the mechanism of tumour immune escape based on the T cell dysfunction and infiltrating in tumour. The TIDE score was used to assess immune escape and immunotherapy of HNSCC patients (http://tide.dfci.harvard.edu/). After that, the CIBERSORT algorithm applied to measure the immune cell abundance of different risk score HNSCC patients in cancer microenvironment with default parameters ^[24]. We evaluated the sensitivity of five common chemotherapeutic agents for high and low-risk patients through pRRophetic package. The immunoassay results are presented in heatmap format by Heatmap package.

2.8.   Statistical analysis

Differences between continuous variables were compared using the student's t-test or the Mann-Whitney test, and differences between categorical variables were compared using the chi-square test using the limma package. COX regression analysis was used to screen for prognostic lncRNAs, while Lasso regression analysis was used to screen for stable prognostic models to avoid overfitting. ROC curves were utilized to evaluate the ability to predict OS of HNSCC patients. Kaplan-Meier analysis was performed to elaborate survival probability in different risk score patients. The independent prognostic signature in risk score and clinical-pathological characteristics was determined using multivariate and univariate COX regression models. Significant differences were considered to exist when P-value < 0.05.

3.   Results

3.1.   Cuproptosis-related gene expression profile

Figure 1 presents the study flow. Initially, the 17 cuproptosis-related genes were gathered from previous studies. The UMAP algorithm descends the multidimensional single-cell sequencing data and projects the 16 cell clusters into a two-dimensional plane (Figure S1A). The 16 cell clusters were manually annotated into seven cell subgroups based on marker genes obtained from previous studies (Figure S1B and Figure S2A). One of these clusters, Figure 2A shows the high expression of 13 cuproptosis-related genes in the cancer cluster the lipoic acid pathway (FDX1, LIAS, LIPT1, DLD), the pyruvate dehydrogenase (PDH) complex related genes (DLAT, PDHA1, PDHB, MTF1, GLS, CDKN2A), protein lipoylation enzymes related genes (DBT, GCSH, DLST, DLAT). In addition, the genes (SL3CA1, ATP7A, ATP7B) encoding the copper transporter protein were highly expressed in endothelial cell cluster, while NLRP3 highly expressed in the macrophage cluster. 12 DEGs were shown to aggregated between tumour and normal samples in the heatmap (Figure 2B). Subsequently, we explored the enrichment functions of the DEGs in three categories, showing that DEGs were enriched in the alpha-amino acid metabolic process, biosynthetic process, mitochondrial matrix and oxidoreductase activity acting on the aldehyde or oxo group of donors (Figure 2C). Subsequently, DEGs were associated with citrate cycle, carbon metabolism, and pyruvate metabolism, according to the KEGG analysis (Figure 2D).

Figure 2.  Cuproptosis-related genes expression profile, enrichment function and pathway analysis. (A) The UMAP plot showed 17 cuproptosis-related gene expression profiles in the microenvironment. The first UMAP plot shows the results of the cell cluster annotation. The colour change from purple to yellow represents an increase in the cell cluster's density of cuproptosis-related gene expression. (B) Heatmap for differentially expressed cuproptosis-related genes. (C) The GO enrichment function analysis. (D) The KEGG pathway analysis.

DownLoad: Full-Size Img PowerPoint

3.2.   Construction and testing of the cuproptosis-related lncRNA prognostic model

According to the annotation file on the Genecod website, we screened 16,876 lncRNAs in bulk sequence data of HNSCC patients and 486 cuproptosis-related lncRNAs are presented in supplementary file (Table S2). Univariate COX analysis confined 33 significant prognostic lncRNAs (Figure 3A). The lowest point of the curve showed 22 as optimal turning parameters through Lasso regression, and its coefficient were the crossover values of the variables and parameters (Figure 3B, C). Multivariate COX regression analysis was screened out 11 lncRNAs constructing lncRNA prognostic model (MIR9-3HG, AC004943.2, AL133395.1, AL139287.1, AC005332.1, RAB11B-AS1, AC108693.2, AL160314.2, AC144652.1, AC005076.1, AC025171, Figure 3D). Coefficient value of prognostic model is shown in Figure 2E, the equation of which is as follows: risk score = -0.299634761231868 * MIR9-3HG + 0.416615078213639 * AC004943.2 + 0.406186176575351 * AL133395.1 + -0.32630458882232 * AL139287.1 + -0.996673062991142 * AC005332.1 + -0.529561147888853 * RAB11B-AS1 + -0.685936586741667 * AC108693.2 + 3.82815712962467 * AL160314.2 + 0.262183908147089 * AC144652.1 + 0.335576038606579 * AC005076.1 + -0.629093407843113 * AC025171.2. The heatmap elaborated the relationship between DEGs and 11 prognostic lncRNAs (Figure 3F). The model was separated into training and test sets, where patients were split into high and a low-risk groups depending on median risk score. Distinct separation between two risk score patients were shown in PCA and TSNE plots (Figure 3G, H). It was clearly observed that the OS rate was worse in the group with higher risk when compared to the group with lower risk in each set (Figure 4A–C). The scatter map demonstrates gradual decrease in survival status of HNSCC patients as the risk score increases in all sets (Figure 4D–F). The AUC for 1-, 3-, and 5-years OS were 0.742, 0.786, 0.722 in training set, 0.658, 0.654, 0.602 in test set, and 0.679, 0.713, 0.656 in whole set (Figure 4G–I).

Figure 3.  Constructing cuproptosis-related lncRNA prognostic model. (A) A forest map plotting 33 cuproptosis-related lncRNAs in univariate COX regression. (B) The optimal turning parameters (log λ) are selected through the thousandfold cross-validation. (C) Lasso coefficients for 33 cuproptosis-related lncRNAs. (D) A forest map plotting 11 cuproptosis-related lncRNAs in multivariate COX regression. (E) Coefficient cuproptosis-related lncRNA prognostic model regression. (F) Correlation of prognostic lncRNA with cuproptosis-related genes. PCA (G) and TSNE (H) demonstrates the distribution of different risk score patients.

DownLoad: Full-Size Img PowerPoint

Figure 4.  Testing of the cuproptosis-related lncRNA prognostic model. (A–C) An overview of Kaplan-Meier curves for training, test, and the whole set, respectively. (D–F) Trends in survival status as risk score increases for training, test, and the whole set, respectively. (G–I) ROC curves of the prognostic model for predicting OS of HNSCC patients in training, test and the whole set, respectively.

DownLoad: Full-Size Img PowerPoint

3.3.   Identification of an independent prognostic signature

We identified valuable prognostic signatures in conjunction with clinical-pathological characteristics for testing the role of prognostic signatures. Table S3 shows the clinical baseline profile for two risk score patients HNSCC patients. The model may serve as an independent prognostic signature compared to other clinical-pathological characteristics through univariate and multivariate COX regression in each set (Figure 5A–C). A C-index and ROC cure were evaluated across the whole set for assessing the diagnostic value and capability of each signature. C-index demonstrated that the predictive capability of the model was better than additional prognostic signatures like age and stage (Figure 5D). It was found that the AUC values of the model at 1-, 3-, and 5-years were 0.679, 0.713, 0.656 respectively and higher than clinical features, thereby testing the prognostic model. (Figure 5E–G). The nomogram combining clinical-pathological characteristics and the model might accurately predict the OS of HNSCC patients (Figure 6A). The relationship between prognostic lncRNA and clinical pathological characteristics were shown in heatmap (Figure 6B).

Figure 5.  Identification of independent prognostic signature. (A–C) Forest plots for risk score and clinical pathological characteristics correlated with prognosis. (D) C-index curve demonstrated the accuracy of the model and clinical pathological characteristics for prognostic prediction in the whole set. (E–G) ROC curves evaluated the prediction ability of model and clinical-pathological characteristics for HNSCC patients in the whole set.

DownLoad: Full-Size Img PowerPoint

Figure 6.  (A) A nomogram for prediction OS of HNSCC patients according to clinical pathological characteristics and the model. (B) Heatmap for 11 cuproptosis-related prognostic lncRNAs and clinical pathological characteristics.

DownLoad: Full-Size Img PowerPoint

3.4.   Enrichment function analysis and immunogenomic landscape analyses

Since cuproptosis-related lncRNA prognostic models have been tested and can be considered as independent prognostic factors. Therefore, we performed a differential analysis of the whole set of patients with HNSCC between different risk groups. Differentially expressed genes between the two risk patients were screened with the FDR value > 1 and cut-off of P-value < 0.05. Subsequently, we performed enrichment function analysis of differential genes and pathway analysis of DEGs by GO and KEGG analysis. The results of GO enrichment demonstrated that DEGs were mainly enriched in immune function and cell structure, such as immune response-activating signal transduction, immune response-activating cell surface receptor signaling pathway, T cell receptor complex, immunoglobulin complex, circulating (Figure 7A). KEGG analysis demonstrated DEGs enriched in the agent metabolism-cytochrome P450 pathway (Figure 7B). The enrichment analysis results suggest that immune responses and infiltration play an essential role in cancer development. The TIDE score was used to signature cancer escape mechanisms, thus, high-risk score patient of HNSCC patients is potentially benefited from immune checkpoint inhibitor therapy (Figure 7C). The effect of immune function and infiltration of HNSCC patients was obtained by scoring immune cell marker genes through multiple algorithms. The ssGSEA indicated the T cell function of low-risk patients was mainly enriched in cytolytic, inflammation-promoting, HLA, T cell stimulation and inhibition (Figure 7D). Figure 5E presented the immune subgroup cells, such as monocytes cells, CD4 memory T cells and gamma-delta T cells infiltrated the cancer microenvironment of low-risk patients.

Figure 7.  Immunoassay for the prognostic model. GO enrichment function (A) and KEGG pathway analysis (B) of DEGs. (C) TIDE score for two risk score patients. Heatmap of immune function (D) and infiltration (E) in the high and low-risk score patients.

DownLoad: Full-Size Img PowerPoint

3.5.   chemotherapeutic agent sensitivity

To explore potential treatment strategies in prognostic models, we explored the sensitivity of five common chemotherapeutic agents (IC50) for treating different risk score groups. The IC50 value of docetaxel and cisplatin for low-risk score patients was significantly higher than high-risk score, which indicates that these chemotherapeutic agents were more effective in high-risk score patients (P < 0.05, Figure 8A, B). Nonetheless, the IC50 for gemcitabine, paclitaxel and methotrexate were not statistically different in the two risk score patients. (Figure 8C–E).

Figure 8.  Five common chemotherapy drug sensitivity. The IC50 value of docetaxel (A) and cisplatin (B) were significant differences between the high and low-risk score patients but not in gemcitabine (C), paclitaxel (D), and methotrexate (E) *P < 0.05; **P < 0.01; ***P < 0.001.

DownLoad: Full-Size Img PowerPoint

4.   Discussion

Cisplatin is a main anti-cancer agent used in the treatment of HNSCC, which acts by causing DNA damage and inducing apoptosis in cancer cells. However, cancer cells develop resistance to cisplatin as treatment progresses, reducing the treatment effectiveness and finally leading to treatment failure because ATPase copper transporting beta (ATP7B) encloses the cisplatin within vesicles in the extracellular space to prevent cisplatin-induced apoptosis ^[25,26]. Cancer cell death is induced by altering copper ion concentrations in agents or by changing APT7B expression, providing a novel insight into the treatment of HNSCC ^[27,28]. In this study, ATP7B is predominantly highly expressed in endothelial cell clusters which structure the tumour microenvironment and promote cancer development. Cuproptosis is copper-induced cell death by targeting lipoylated tricarboxylic acid cycle proteins ^[6]. KEGG analysis identified the TCA cycle as a potential pathway of HNSCC. The pyruvate dehydrogenase (PDH) complex controls pyruvate entry into the TCA cycle and inhibition of the PDH complex could promote the Warburg metabolic and malignant phenotype in HNSCC through inhibiting phosphorylation of PDH ^[29]. The PDH complex in cuproptosis includes DLAT, PDHA1, PDHB, MTF1, GLS, and CDKN2A, which are highly expressed in the HNSCC cancer cluster, which demonstrate a novel perspective for controlling the development of cancer. In addition, the lipoic acid pathway, a key mediator of copper-induced cell death, is highly expressed in clusters of HNSCC cells.

Genome-wide studies have demonstrated that alterations and mutations of a substantial number of lncRNAs may promote cancer occurrence and metastasis ^[30]. We explored 33 cuproptosis-related prognostic lncRNAs from TCGA dataset and constructed a cuproptosis-related lncRNAs prognostic model based on 11 prognostic lncRNAs (MIR9-3HG, RAB11B-AS1, AC004943.2, AL133395.1, AL139287.1, AC005332.1, AC108693.2, AL160314.2, AC144652.1, AC005076.1, AC025171.2). miRNA could positively regulate the process of epithelial-mesenchymal transition in HNSCC cells, with miR-9 maintaining HNSCC cell invasion and growth in response to radiotherapy-induced cell death ^[31,32]. Therefore, miR-9 is a biomarker of radiotherapy efficacy to guide physicians in the optimal treatment plan for HNSCC patients. MIR9-3HG, encapsulated miR-9 gene, was up-regulated in HNSCC tissue compared to control group and associated with the prognosis ^[33]. Shao Lina found that RAB11B-AS1 was down-regulated in HNSCC cancer tissue and is an independent prognostic signature of HNSCC ^[34]. RAB11B AS1 promoted the progression of hepatocellular carcinoma by downregulating in m6A dependent manner via METTL16 ^[35]. However, the regulatory mechanisms and functions of other lncRNAs in HNSCC are still unknown, and more research is needed to clarify them. In this study, cuproptosis-related lncRNAs are independent risk factors for HNSCC prognosis. The lncRNAs prognostic model predictive performance and diagnostic value for 1-, 3-, and 5-year OS is significantly higher than other clinical characteristics. To make it possible to apply the predictive model more easily, a nomogram was constructed to score patients to evaluate their 1-, 3-, and 5-year OS. Although many studies in the mechanism of action of lncRNAs in HNSCC are underway, prognostic studies of lncRNAs in HNSCC patients are still essential.

Survival status for patients with high-risk score were worse than those with low-risk score. It is therefore essential to select the correct treatment and assessment strategy for cancer according to the different risk score patients. In this study, GO enrichment functional analysis may point us to potential therapeutic strategies for HNSCC patients. GO enrichment analysis of DEGs revealed molecular functional enriched in immune response-activating signal transduction and immune response-activating cell surface receptor signaling pathway. This is strong evidence that treatment of patients in different risk score groups require improvement of the immune environment or pathways. Both traditional chemotherapy and immunotherapy for HNSCC involve modulation of the immune system and activation and suppression of immune cells such as cytotoxic T cells and helper T cell ^[36]. Chemotherapy is generally considered to have immunosuppressive effects due to myelosuppressive effects, resulting in lymphopenia and neutropenia. However, traditional chemotherapeutic drugs such as cisplatin can affect the immune system to varying degrees and thus enhance the anti-tumor immune response ^[37]. PD-1 is a co-inhibitory receptor on the cell surface of cytotoxic T lymphocytes. Ligation of PD-1 and PD-L1 on tumor cells or antigen presenting cells (APCs) suppresses the immune response, which reduced effector and memory T cell production and increased Treg and exhausted T cell production ^[38]. Therefore, blocking this ligation becomes a theoretical basis for anti-PD-1 immunotherapy for the treatment of HNSCC. Recently, PD-1/L1 inhibitors have become a new research direction in immunotherapy and have shown promising efficacy and controlled safety in several oncology trials ^[39]. This is consistent with our results of immunogenomic landscape analysis for different risk score patients. We analyzed immune function and immune infiltration for HNSCC patients and found that activation and suppression of different subtypes of T cells correlated significantly with the patients of different risk groups. Thus, altered survival in HNSCC patients may benefit from altered T cell subpopulation cell infiltration and activation of toxic T cells. Traditional cytotoxic chemotherapy may become increasingly important as an adjunct to immunotherapy for patients with tumours poorly infiltrated with lymphocytes. Esophageal squamous cell carcinoma mice treated with cytotoxic chemotherapy showed significantly elevated amount of CD4 and CD8 T cells in the tumour microenvironment ^[40]. Balermpas et al. demonstrated the prognosis of patients treated with radiotherapy at different sites, where patients with well-infiltrated CD8 and CD3 tumour-infiltrating lymphocytes had a more favourable prognosis than those with poor infiltration, with no differences in other clinicopathological factors ^[41]. This is also consistent with our results of significantly lower T lymphocyte infiltration and OS in the high-risk HNSCC patients. In the present researches, the IC50 of high-risk score patients were significantly lower in docetaxel and cisplatin than in the low-risk score patients. This reveals that the high-risk score patients are more sensitive to docetaxel and cisplatin chemotherapeutic agents. More studies are needed to determine whether the treatment effect can be improved by altering the T-cell infiltration in patients at a later stage. Many researches have indicated that grade of immune infiltration is related to survival status of patients with high CD8T cell ratios showing better survival outcomes in both HPV-positive and -negative patients ^[40,42]. Immune checkpoint blockade could specifically reverse the suppression of tumour-infiltrating CD8+ T cells and promote natural killer cells immune function through blocking PD-1/L1 receptors ^[43,44]. Although PD-1/L1 blockers have shown significant success in many tumours, not all patients benefit from them. Yajuan Zhang found that ROS clearance in the TCA cycle and NF-κB activation collectively promote tumour cell survival in response to inflammatory factor stimulation, enhance tumour cell tolerance to cytotoxic T lymphocytes and ultimately promote tumour immune escape ^[45]. In our study, high-risk patients with lower TIDE scores are less likely to experience immune escape and could benefit from immunotherapy. However, altering ROS levels in the TCA cycle to promote copper-induced cancer cell death and combined immunotherapy still requires more experimental studies.

Although we have constructed a stable prognostic model and predicted target immunotherapy loci, our study has several limitations. First, we tried to find additional independent HNSCC sequencing data in other public databases but only the TCGA database provided complete clinical-pathological characteristics and bulk RNA sequencing data. Secondly, we need to expand the clinical case data and iteratively test and modify the prognostic model. Thirdly, researchers need to discover the regulatory mechanism governing the function of lncRNAs and cuproptosis-related genes in HNSCC. However, I believe that with the following collection of HNSCC patients at our hospital and the discovery of more theoretical evidence to support this, we will work forward based on this study.

In summary, this study recognizes that cuproptosis-related genes were closely associated with HNSCC cell clusters in the tumour microenvironment. A cuproptosis-related lncRNAs prognostic model combined clinical-pathological characteristics was constructed and used to explore the immune landscape and tumour agent resistance in HNSCC patients, providing new comprehensive loci for immunotherapy and antitumor agent development.

Conflict of interest

The authors declare there is no conflict of interest.