Regulatory Role of IL-6 in Immune-Related Adverse Events During Checkpoint Inhibitor Treatment in Melanoma

1. Introduction

Metastatic melanoma is an advanced and aggressive form of skin cancer that arises from the uncontrolled growth of pigment-producing cells called melanocytes [1]. Immune checkpoint inhibitors (ICI) therapies have emerged as a revolutionary approach in the treatment of melanoma. ICI therapies focus on modulating the immune system to enhance its ability to recognize and attack cancer cells. Various ICI has been discovered to inhibit specific proteins (e.g., cytotoxic T-lymphocyte-associated antigen 4 (CTLA-4) [2], programmed cell death protein 1 (PD-1) [3] in tumour immune environment and thereby eliminate of tumour cells [3,4,5,6]. While ICI therapies hold great promise in cancer treatment, they present a distinct challenge by potentially inducing autoimmune phenotypes [7]. The mechanisms that enable the immune system to target cancer cells are also supposed to inadvertently lead to immune-related adverse events (irAEs) responsible for the exacerbation of autoimmune diseases [8]. This dual nature of immune checkpoint therapies underscores the delicate balance that must be struck between activating the immune system to combat melanoma and preventing it from attacking the body's tissues resulting in irAEs [9,10,11,12,13,14]. Understanding the cross-talk between melanoma and autoimmune diseases is challenging due to the involvement of a large number of immune cells and is vital for the safe use of ICI therapies in the regulation of tumor growth. In this manuscript, we approached this challenge using network-based approaches by constructing and analysis molecular and signaling network associated with the tumor-immune cross-talk. Detailed work flow of research is given in Figure 1.

2. Results and Discussion

2.1. Protein-Protein Interaction Network at the Interface of Melanoma and Autoimmune Diseases

We extracted genes associated with MM (n=504), RA (n=2722), UC (n=1458), and CD (n=1382) from the DisGenNet database. We found a total of 132 common genes in all four disease phenotypes for which a protein-protein interaction (PPI) network was prepared using the STRING database (Figure 2). We considered these common genes as the connecting links between melanoma and investigated autoimmune diseases.

Out of 132 genes used for the construction of the PPI network, 25 genes did not exhibit any connections with other genes, and hence a network among the common 107 genes (Figure 3).

2.2. Identification of Hub Genes from the PPI Network Associated with the Crosstalk between Melanoma and Autoimmune Diseases

To identify key hub genes in the PPI network prepared using the common genes in melanoma and autoimmune diseases, we used the Cytoscape plugin Molecular Complex Detection (MCODE). The MCODE algorithm detects interconnected network clusters based on a k-core score that represents the maximal number of connected subgraphs with all the nodes connected by a minimum number of degrees k [15]. The MCODE algorithm detected six highly connected subnetworks represented as MCODE modules [Table 1]. For the identification of hub genes responsible for the crosstalk between melanoma and autoimmune diseases, we selected module 1, containing 16 genes with the highest MCODE score of 10.13.

2.3. Pathway Enrichment Analysis of top MCODE Cluster

To identify biological processes and pathways that might play a key role in the crosstalk of melanoma and autoimmune diseases, we performed pathway enrichment analysis of genes associated with the top MCODE module using the Reactome database 2022 in Enrichr web-based server (https://maayanlab.cloud/Enrichr/). All the enriched pathways (Figure 4), along with the pathway p-values, adjusted p-values, and gene/protein sets for each case (Supplementary Table S1).

Interleukin-4 (IL-4) and Interleukin-13 (IL-13) signaling pathways, which play key roles immune regulation and inflammation and are primarily associated with allergic responses and Th2 immunity [16], were among the top enriched pathways. Previous studies also suggest that these pathways shape the tumor microenvironment and promote tumor progression by modulating immune responses and survival pathways [17]. In autoimmune diseases like RA and UC, IL-4 and IL-13 contribute to excessive inflammation and tissue damage, fostering autoantibody production and B cell survival[18]. Besides, we also found Interleukin-6 (IL-6), Interleukin-10 (IL-10), and Interleukin-1 (IL-1) signaling cascades among the top enriched pathways that regulate melanoma [19,20] and autoimmune diseases [21], exerting diverse effects on inflammation, immune dysregulation, and disease progression. Elevated IL-6 levels in melanoma correlate with advanced disease stages and therapy resistance, while in autoimmune diseases, IL-6 drives inflammation and tissue damage [21]. Conversely, IL-10 exhibits dual roles, suppressing anti-tumor immunity in melanoma yet mitigating inflammation in autoimmune diseases. In melanoma, IL-1 orchestrates tumor growth, angiogenesis, and metastasis [22] by instigating pro-inflammatory cytokine production, fostering melanoma cell invasiveness, and modulating the tumor microenvironment's immune cell composition that may contribute to autoimmune disease pathogenesis, fueling inflammation and tissue damage in conditions like RA [23] and IBD [24] by inciting cytokine production and immune cell activation. Additionally, we also found matrix metalloproteinases (MMPs), collagen degradation, MAPK signaling, and CD163-mediated anti-inflammatory responses among the top enriched pathways that are known for regulating melanoma, and autoimmune diseases [25,26,27,28,29,30,31]. In melanoma, collagen degradation promotes tumor invasion and metastasis by MMP-mediated breakdown of the extracellular matrix (ECM), facilitating melanoma cell infiltration into surrounding tissues and distant metastasis [49,50]. Elevated levels of collagen degradation products also instigate immune-mediated tissue damage and inflammation in autoimmune diseases like RA [34]. The CD163-mediated anti-inflammatory responses potentially affect the inflammation resolution phase resulting in the progression towards chronic inflammatory phenotypes together with the excurbation of autoimmune symptoms. Understanding the intricate interplay of these signaling pathwys is crucial for developing targeted therapies that effectively modulate immune responses together with management of tumor by immune checkpoint inhibitors.

2.4. Identification of Lead Molecule and Molecular Docking

IL6 (Figure 4) was found to regulate five pathways among eight enriched pathways associated with MM and autoimmune disease. This directed network suggests that inhibition of IL6 will reduce the activity of pathways such as IL6, IL4 / IL13 signaling, collagen degradation, CD163-mediating responses, MAPK signalling, and activate IL-10 signaling. Thus, inhibition of IL-6 will not only be able to suppress melanoma metastasis but reduce autoimmune phenotypes, which may exacerbate during immune checkpoint therapies [35,36,37].

To identify potential inhibitors for the IL6 protein, we utilized the FDA-approved library using the Zinc15 Database (https://zinc15.docking.org/), which contained 1615 compounds and the natural compound library of the COCONUT database (https://coconut.naturalproducts.net/), which contained 407,270 unique compounds [38]. We filtered the COCONUT database libraries prior to the virtual screening with IL6 for Lipinski's rule of five[39] in order to consider only drug-like molecules. Only 272,001 compounds were able to pass Lipinski's rule of five filtering criteria.

The active site of IL-6 was selected considering amino acid residues Phe74, Phe78, Leu178, Arg179, and Arg182, which play a significant role (hotspot residues) in its interaction with IL6R [40]. The ‘Libdock’ tool in DS2022, which is a rigid-based docking program that calculates hotspots for the protein using a grid placed into the binding site polar and apolar probes [41,42], was used for the screening process of drug libraries. The top 20 compounds based on the libdock score were further considered for flexible docking analysis using the ‘CDOCKER’ tool in the DS2022 [Table 2]. In Flexible CDOCKER docking, the compound is allowed to adapt its conformation to fit the binding site of the receptor better, and the receptor can also undergo conformational changes to accommodate the compound. This flexibility is particularly important in accurately predicting the binding affinity and mode of interaction between a compound candidate and its target protein. Only nine compounds out of 20 could be further docked with IL6 using CDOCKER protocol.

We selected the top two compounds (Figure 5) because these have the highest CDOCKER energy value (CNP0003038: -41.6684 kcal/mol) and (ZINC03809192: -34.7136 kcal/mol), and one belongs to the natural compound library and second from the ZINC database.

We performed the literature survey for these two compounds for their toxicity, bioassays, and role in tumor / autoimmune disease regulation. We found that the compound CNP0003038, also referred to as 1,1'-ethylenebis-L-tryptophan (EBT; PubChem CID: 3905118), enhances the proliferation of EoL-3 eosinophilic leukemia cells, induces the release of eosinophil cationic protein from isolated human peripheral blood eosinophils resulting in eosinophilia-myalgia syndrome [43]. EBT is also shown to induce IL-5 production in isolated human T cells and induces inflammation, mast cell infiltration, fascia thickening, and adipose tissue fibrosis in an eosinophilia-myalgia-syndrome mouse model [44]. Due to the possible toxic effects of compound CNP0003038 on the immune system, we removed it from our further analysis.

The compound from the ZINC database ‘ZINC000003809192’, also known as amprenavir, is primarily known as a protease inhibitor and used in the treatment of HIV/AIDS. It inhibits the HIV protease enzyme, thereby blocking the cleavage of viral polyproteins into functional proteins, ultimately hindering viral replication [45,46]. Amprenavir was also included in the investigation of FDA-approved small molecule drugs through in-silico screening, assessing their potential as inhibitors of extracellular signal-regulated kinase (ERK) and apoptosis inducers in MCF-7 human breast cancer cells [47]. Based on all the above facts, we selected only Amprenavir for the molecular dynamic simulation.

We employed molecular docking and molecular dynamics simulation analyses to investigate the interaction patterns, stability, and flexibility of the docked complex, which helped us explore the interaction and stability of Amprenavir with IL-6 during the simulation. Our molecular docking studies suggest that Amprenavir forms three hydrogen bonds and four hydrophobic bonds with the IL6 amino acid residues SER176, CYS73, MET67 and, ARG179, LYS54, LYS66 (more information on the bonds is provided in the supplementary Table S2). To further check if Amprenavir interferes with the binding of IL6 with IL6R, we performed additional protein-protein docking using the HDOCK tool [48]. For this, we first used IL6 and IL6R complex (PDB ID: 1P9M) [49] and redocked the protein units using HDOCK as a control scenario. Next, we used the IL6 in complex with ‘Amprenavir’ and performed the protein-protein docking with IL6R again using the HDOCK tool. Both these scenarios were compared with each other to evaluate the effect of ‘Amprenavir’ on IL6 and IL6R interactions. We observed that the docked complex of IL6_IL6R has a higher docking energy of -398.1 kcal/mol thus more stable in comparison to the IL6_Amprenavir_IL6R complex (-262.02 kcal/mol). We further compared the impact of ‘Amprenavir’ on the number of bonds formed between IL6 and IL6R.

The residues essential for IL6 binding to IL6R, includes IL6: LYS27, GLN28, ARG30, PHE74, PHE75, PHE78, LEU178, ARG179, ALA180 and ARG182 [49] which form various bonds with IL6RA: GLU163, GLN190, PHE229, ASP253, GLU277, GLU278, PHE279, and IL6RB: LYS118, LYS119, ARG128, VAL167, TYR168, PHE169, VAL230, VAL264 residues (more information on the bonds formed between IL6 and IL6R is provided in the supplementary Table S3).

Our analysis revealed that the IL6_ Amprenavir_ IL6R complex lost 6 significant bonds that were formed between the IL6_IL6R complex. However, one new bond formed between IL6 and IL6R in the presence of amprenavir (IL6:ARG179 with IL6RA: GLU163) (supplementary Table S3). We also observed IL6 amino acid residues LYS66, SER176, ARG179, and IL6R residues GLY164 and CYS192 forming bonds with amprenavir (Figure 6 a and b) (supplementary Table S3). The energies and bond assessments of the docked complex show that the compound amprenavir functions as an inhibitor for IL6 interaction with IL6R, disrupting numerous bonds originally formed between the protein and its receptor.

2.5. Molecular Dynamics Simulation

To evaluate the IL6_amprenavir docked complex flexibility and overall stability through time-dependent molecular dynamics (MD) simulation, we used the ‘Standard Dynamics Cascade’ protocol of DS 2022. The analysis of the simulation involved the use of RMSD (Root Mean Square Deviation), RMSF (Root Mean Square Fluctuation) and Radius of gyration (Rg) (Figure 7). We observed the RMSD of the docked complex (IL6_Amprenavir) for 50 ns; and found that the complex achieved convergence and stabilizes at around 20 ns (Figure 7 b).

We comprehensively analysed residue fluctuations in RMSF, which are crucial for IL6 binding with IL6R. Residues within the hydrophilic domain (Lys27, Arg30, Phe78, Arg179, and Arg182) form salt bridges with IL6R proteins. In the molecular docking of the IL6_Amprenavir_IL6R complex, we observed the disruptions of salt bridges involving residues Lys27 and Arg182 between IL6 and IL6R. The same happens in Phe 78, disengaged from its bond with IL6R (supplementary Tables S3). All these residues take part in the binding of the amprenavir into the cavity of IL6. During the MD simulation, all these residues showed minimal fluctuations of 0.4 (Å), which indicates that these IL6 residues are tightly engaged with the drug ‘amprenavir’ and are not available for interaction with IL6R. Furthermore, no residues showed a fluctuation of more than 0.7 Å (Figure 7c). In the case of the radius of gyration, during the MD run after the binding of the Amprenavir, IL6 started to achieve a more compact structure (Figure 7d). The comparison of bonds formed between IL6_amprenavir complex before and after the MD simulation is shown in Figure 7a and in the supplementary Tables S2 and S4. The analysis of the final pose of IL6_Amprenavir complex after MD simulation indicates an increase in bonds compared to the initial docked pose. Initially, it formed three Hydrogen bonds with various IL6 residues which increased to five after a 50 ns production run. Overall, our analysis indicates that the compound amprenavir binds with IL6 and forms a stable complex at the same site that is associated with IL6R interaction.

3. Methods and Methodology

3.1. Data Collection and Protein-Protein Interaction (PPI)

The genes related to MM, RA, UC & CD were extracted from the DisGenNet database (https://www.disgenet.org/), which is one of the largest publicly available collections of genes and variants associated with human diseases [50]. The data extracted from DisGeNet was subsequently analysed to identify common genes among the above-mentioned disease phenotypes. All the common genes were further used to explore their protein-protein interactions using the STRING database (http://string-db.org) [51] with a confidence score cutoff of 0.70.

3.2. Identification of Highly Interconnected Clusters in the Tumor-Autoimmune PPI Network

We used the Cytoscape plugin Molecular Complex Detection (MCODE) [53,54,55] to identify highly connected clusters in the protein-protein interaction network of the common genes associated with tumour and selected autoimmune disease phenotypes. The degree cutoff = 2 was used to control the genes to become part of the cluster. New members are added to the cluster only if their node score deviates from the cluster’s seed nodes score by less than the set cutoff of 0.2. The k-score cutoff = 2 was used to filter out clusters that do not contain a maximally inter-connected sub-cluster of at least a degree score of 2. The max. depth of 100 was set to limit the distance from the seed node within which MCODE can search for cluster members. This approach allowed us to partition the network based on its topology, pinpointing densely connected regions within the protein-protein interaction network associated with the crosstalk between melanoma and immune checkpoint therapies inducing irAEs.

3.3. Pathway Enrichment Analysis

For the top cluster identified using MCODE analysis, we identified enriched pathways using Enrichr (https://maayanlab.cloud/Enrichr/) [52]. We specifically focused on the pathways that are present in the Reactome 2022 database with the significance threshold of p-value (0.05) for the enrichment analysis. The pathways were sorted on the basics of the combined score [53], which is described as

c = log(p)*z,

where c is the combined score, p is the Fisher exact test p-value, and z is the z-score for deviation from expected rank. Further, we filtered enriched pathways that were previously identified to play a role in melanoma and autoimmune diseases. From the filtered pathways, we identified nodes that are common in all the enriched pathways for future analysis.

3.4. 3D Structure Preparation and Screening of Lead Compounds

We conducted virtual screening and molecular docking analyses to elucidate the inhibition mechanism and identify potential lead compounds for IL6 inhibition. To perform these analyses, we use a 3D model of IL6 (PDB ID: 1ALU) from the RCSB PDB database (https://www.rcsb.org/) [54]. The model was subjected to the ‘Prepare Protein’ protocols of the Biovia Discovery Studio 2022 (DS2022) [55] using the CHARMm force field [56]. The structure was further optimized using the ‘Smart Minimiser’ algorithm to achieve a stable state through energy minimization. The minimization process was done in 2000 steps, with an RMS gradient tolerance set at 0.1. With the help of the ‘Receptor-ligand Interaction’ tool of DS2022, we defined the IL6 binding site based on the amino acid residues that play a key role in its interaction with IL6R [54]. To identify potential inhibitors for IL6 proteins, we conducted a virtual screening of natural products using the Collection of Open Natural Products (COCONUT), an extensive and well-annotated resource for natural products [38], Additionally, we included an FDA-approved drug library from the ZINC 15 database[57] for the virtual screening using the ‘LibDock’ protocol of the DS2022. All the docked compounds were further subjected to the Flexible ‘CDOCKER’ program in the DS2022 [55,58]. Finally, we performed the protein-protein docking between IL6R and IL6 in complex with the screened drug molecules using the HDOCK tool [48].

3.5. Molecular Dynamic Simulation

We assessed the binding affinity of the IL6 inhibitor docked at the IL6R interaction site, we performed MD simulation using DS2022. The MD simulation was performed in an implicit solvent environment to investigate the stability, conformational changes, and dynamic behavior of the inhibitor in the binding cavity by examining the formation of diverse electrostatic interactions. For the MD simulation, the complex was subjected to an initial minimization phase consisting of 1000 steps using the steepest descent algorithm, followed by an additional 2000 steps employing the conjugate gradient method with the CHARMM force field [56]. After minimization, the systems underwent a heating phase where the initial temperature was incrementally increased from 50 K to 300 K in 50 ps intervals. Subsequently, an equilibration step lasting 100 ps was performed. The adjusted velocity frequency was configured at 50 for the heating and equilibration phases. Subsequently, a production run of 50 ns was conducted within an NVT assembly (maintaining normal volume and temperature) at a constant temperature of 300 K, with results being saved at intervals of 0.02 ns. For the entire simulation run, we analyzed trajectories consisting of 25,000 conformations. Various properties, including root-mean-square deviation (RMSD), root-mean-square fluctuation (RMSF), and the radius of gyration (ROG), were examined using the 'Analyze Trajectory Protocol' of DS2022

4. Conclusions

The introduction of ICI, such as PD-1, PD-L1, and CTLA-4 inhibitors, has revolutionized the clinical management of metastatic melanoma (MM). These inhibitors have shown remarkable efficacy in controlling the growth of tumors. However, their use also led to irAEs, resulting in exacerbating autoimmune diseases such as RA, UC, and CD in melanoma patients. Our research has explored the interface between MM and autoimmune diseases, aiming to identify potential druggable targets and understand the cross-talk between these conditions for the safe use of ICI in MM management. Using an integrative network approach, IL6 was identified as a promising target at the interface between MM and autoimmune diseases. Through structure biology approaches, lead compounds capable of inhibiting IL6, such as Amprenavir, have been identified. However, further laboratory experiments are needed to validate the efficacy of Amprenavir together with ICI.