Cocaine addiction-related data GSE54839 (Homo sapiens), GSE67281 (Homo sapiens), GSE186981 (Mus musculus), and GSE155313 (Mus musculus) were downloaded from the Gene Expression Omnibus (GEO) database (http://www.ncbi.nlm.nih.gov/geo). The dataset GSE54839 was chosen as our experimental set, which is based on the GPL6947 platform (Illumina HumanHT-12 V3.0 expression beadchip). This microarray-based study determined the profiles of midbrain gene expression in chronic cocaine abusers (n = 10) and well-matched drug-free control subjects (n = 10). Array-related procedures were performed in triplicate for each subject.

GSE67281 (Homo sapiens), GSE186981 (Mus musculus), and GSE155313 (Mus musculus) were chosen as our validation sets. GSE67281 is an expression profile in postmortem human midbrain specimens from chronic cocaine abusers (n = 11) and well-matched control subjects (n = 11). The GSE186981 is RNA-Seq data in hybrid mouse diversity panel (HMDP) mouse strains of nucleus accumbens (NAc) and prefrontal cortex brain regions (PFC). GSE155313 is the RNA-Seq data from the VTA region of mouse that underwent one of four commonly used paradigms: acute home cage injections of cocaine, chronic home cage injections of cocaine, cocaine-conditioning, or intravenous-self administration of cocaine.

Human PPI data were downloaded from the STRING database (https://string-db.org/) (Szklarczyk et al., 2021). With a combined score of >900 as the threshold, a total of 230,524 interactions between 11,763 genes were obtained.

The dataset GSE54839 was differentially analyzed using the R package “limma” to obtain their differential genes. The p-value of <0.05 and |log₂F>p20mm| > 0.2630344 (i.e., fold change ≥ 1.2 or fold change ≤ 0.8) were considered statistically significant.

To obtain the interaction relationships between DEGs, the downloaded PPI data were filtered using the DEGs obtained from the microarray data GSE54839. By using the gene interactions as edges and DEGs as nodes, a differential gene network was constructed. After removing scatters from the network, the core network was defined as a PPI network related to cocaine addiction.

We proposed a centrality algorithm integration strategy to analyze genes in cocaine addiction-related PPI networks. The scores of the node in network under each centrality algorithm were calculated separately by applying a series of centrality measures, including degree, edge-percolated component (EPC), Laplacian centrality, maximum neighborhood component (MNC), Katz radiality, and semi-local centrality (SLC). The intersection of the top 10 genes of each centrality algorithm was considered the key gene.

In this study, G is the cocaine addiction-related PPI network we built, and V(G) is the collection of nodes in the network. For node x in G, N(x) is the set of direct neighbors of x in G. For collection A, |A| is used to represent the number of elements in the collection. The specific algorithms are as follows:

The Human Protein Atlas (HPA; https://www.proteinatlas.org/) database creates a brain-centric knowledge resource on RNA and protein expression in three mammalian brains: human, pig, and mouse (Sjöstedt et al., 2020). The RNA expression of key genes in different brain regions in humans and mice was searched in the brain section of the HPA database.

The Comparative Toxicogenomics Database (CTD, https://ctdbase.org/) was used to obtain associations between key genes and cocaine addiction. In the CTD database, the inference score reflects the degree of similarity between the CTD chemical–gene–disease network and a similar scale-free random network (Davis et al., 2023). The higher the score, the higher the degree of association between the disease and the gene.

The role of key genes in the mechanism of cocaine addiction was identified by text mining in the PubMed database. The search keywords were “cocaine addiction” and the four key genes.

Three independent sets GSE67281 (human), GSE186981 (mouse), and GSE155313 (mouse) were used for pre-addiction and post-addiction differential expression analyses to verify key genes, and the threshold and differential analysis methods were consistent with those of the experimental set GSE54839.

The R package “homologene” was used to search for homologous genes between the human and the mouse. The “homologene” package is a package based on the NCBI HomoloGene (https://www.ncbi.nlm.nih.gov/homologene/) database. The HomoloGene database is a system that can automatically detect congeners in human and mouse genes (NCBI Resource Coordinators, 2014).

To investigate the possible molecular mechanisms of key genes for cocaine addiction, we used the Kyoto Encyclopedia of Genes and Genomes (KEGG; http://www.kegg.jp/ or http://www.genome.jp/kegg/) database for enrichment analysis. An adjusted p-value of <0.05 was considered to be statistically significant.

In order to identify potential therapeutic drugs for cocaine addiction, a network of targeted relationships between nervous system drugs and key genes was constructed. Nervous system drug information was retrieved from the ATC classification system of the Drugbank (www.drugbank.ca) database (Wishart et al., 2018). The targeted effects of key genes with nervous system drugs were reflected in the CTD, where drugs that affect the expression level of genes were our screening criteria. According to the targeting relationship between key genes and nervous system drugs, a targeted relationship network between nervous system drugs and key genes was constructed. Cytoscape was used to visualize this network. Finally, the network was analyzed to screen for potential therapeutic drugs for cocaine addiction.

A total of 724 DEGs were identified from the GSE54839 dataset, including 409 up-regulated genes and 315 down-regulated genes (Supplementary Table 1). The volcano plot was plotted with the “ggplot2” package in R software to visualize the identified DEGs (Figure 2A).

The differential genes were mapped to the downloaded PPI data to obtain a differential gene network consisting of 236 nodes and 316 edges. Removing scatter points in the network, the core network had a total of 153 nodes and 263 edges, which was defined as a cocaine addiction-related PPI network. Cytoscape software was used to visualize the network (Figure 2B).

The scores of each node in the cocaine addiction related PPI network were calculated separately using seven different centrality algorithms (Figure 3A, Supplementary Table 2). Based on our proposed centrality algorithm integration strategy, the top ten genes scored by each algorithm were selected, and their intersections (FOS, IL6, EGR1, and JUN) were regarded as key genes (Figure 3B). The scores of the seven centrality algorithms for the four key genes are shown in Figure 3C.

We performed correlation analyses for four key genes in non-drug control participants (Figure 3D) and chronic cocaine abusers (Figure 3E), respectively. The correlation heat map showed that all four key genes were positively correlated, and the correlation showed a significant increase in the cocaine group. It suggested that these four key genes might be more closely related to each other and had synergistic effects during addiction. In addition, the expression of the IL6 gene was the lowest of the four key genes.

To investigate the expression of four key genes in the brain, we searched the HPA database for the expression of four key genes in different brain regions in the human and mouse (Figures 4A, B). The results showed that four key genes were expressed in all regions of the human brain, and IL6 was expressed in the human midbrain region lower than the other three genes, which is consistent with our findings. IL6 was not detected in the mouse midbrain, hypothalamus, pituitary gland, retina, pons, and medulla.

The correlation scores between each gene in the cocaine addiction-related PPI network and cocaine addiction were searched in CTD. The scores for the top 100 genes are shown in Figure 5. It showed a higher degree of association between the four key genes and cocaine addiction. CTD showed that FOS and EGR1 could be biomarkers of cocaine addiction or play a role in addiction, and EGR1 could be a gene for a therapeutic target in the treatment of cocaine addiction.

The PubMed database showed that four key genes were all associated with cocaine addiction. Both FOS (Fos proto-oncogene) and JUN (Jun proto-oncogene) are members of the AP-1 transcription factor complex. Multiple studies have shown that cocaine affected the expression of FOS proteins (Todtenkopf et al., 2002; Imam et al., 2005; Larson et al., 2010; Lobo et al., 2010). Zhang et al. (2004) prepared CPu extracts from D1 and D3 receptor mutant mice and wild-type control littermates at different time points after cocaine injection and found that ERK activation mediates acute cocaine-induced expression of c-fos (Fos). The study by Xu (2008) found that c-fos might mediate cocaine-induced persistent changes by regulating the formation of AP-1 transcriptional complexes and gene expression. Previous studies have demonstrated that cocaine causes increased expression of the JUN protein (Malaplate-Armand et al., 2005; Paletzki et al., 2008). Cocaine affects the expression of the JUN protein (Imam et al., 2005).

There are some studies proving that EGR1 (early growth active protein 1) and c-fos expressions are reduced after cocaine induction (Helton et al., 1993; Ennulat et al., 1994). In experiments on mutant mice by Valjent et al. (2006), EGR1 was found to play a vital role in cocaine-related behavior. Humblot et al. (1998) found that acute cocaine administration was effective in inducing c-FOS and EGR-1 direct early genes, and cocaine-induced EGR-1 and c-FOS expression was significantly reduced in brain regions of rats.

IL6 (interleukin-6) is a pro-inflammatory cytokine. The study by Halpern et al. (2003) showed that men and women respond weakly to pro-inflammatory challenges to IL6 after intravenous cocaine. In experiments measuring changes in IL6 levels in crack cocaine-dependent adolescents after 21 days of withdrawal, it was found that IL6 was elevated in patients on admission compared to the control group (Pianca et al., 2017).

A differential expression analysis was performed on the human dataset GSE67281, which is the expression profile of human cocaine abusers in the midbrain region. A total of 200 DEGs were identified after annotation, including 110 upregulated genes and 90 downregulated genes. There were 20 intersecting genes in the datasets GSE54839 and GSE67281 (Figure 6), including the key genes JUN, FOS, and EGR1, all of which were downregulated in the addictive state (Table 1), while the IL6 gene was not annotated in GSE67281.

A differential expression analysis and a homology analysis were performed on the mouse dataset GSE155313 from the VTA brain region, and the Fos gene was identified as a downregulated differential gene under four different conditions. The degree of difference in the Fos gene was not the same between chronic and acute home cage injections of cocaine, and even greater in chronic home cage injections of cocaine condition (Table 2). It suggested that the Fos gene plays a crucial role in long-term addiction.

The JUN, FOS, and EGR1 genes were shown to be downregulated differential genes in the human validation set, and the FOS gene was also downregulated in the mouse validation set, which is the same as the experimental set. Analysis of the independent sets showed that the key genes we identified were well-confirmed.

To reveal the roles of the key genes, we performed a KEGG enrichment analysis of all the genes in the cocaine addiction-related PPI network. A total of 75 KEGG pathways were enriched (Figure 7A). The four key genes were mainly enriched in the TNF signaling pathway, cocaine addiction, amphetamine addiction, IL-17 signaling pathway, MAPK signaling pathway, and Toll-like receptor signaling pathway. Previous studies have shown that these pathways were all linked to cocaine addiction (Northcutt et al., 2015; Lewitus et al., 2016; Brown et al., 2018; Ganguly et al., 2019; Montesinos et al., 2020; Bingor et al., 2021). The subnetwork associated with the key genes and their enrichment pathways (Figure 7B) showed that the four key genes were closely connected in the network, among which JUN, FOS, and IL6 were enriched into multiple pathways, and EGR1 was closely related to these pathways and played a very important synergy.

The cocaine addiction pathway and the amphetamine addiction pathway are two enriched addiction-related pathways (Figures 8A, B). They have similar addiction mechanisms, and both have enhanced firing activity of dopamine neurons in the VTA of the midbrain, resulting in enhanced dopamine release from the NAc.

With the stimulation of addictive drugs, the FOS gene induces and maintains an addictive state in the short term of addiction. In human who achieve long-term addiction after further drug use, the FOS gene causes long-term adaptive changes in the brain, and the JUN gene dimerizing with ΔFosB leads to an increased cocaine response. ΔFosB desensitizes c-fos mRNA induction after chronic amphetamine exposure (Renthal et al., 2008). Zhang et al. (2006) have shown that FOS might mediate cocaine-induced persistent changes by regulating AP-1 transcriptional complexes and target gene expression. To sum up, both our key genes JUN and FOS played important roles in the addiction pathways.

To find potential therapeutic drugs for cocaine addiction, a network of targeted relationships between nervous system drugs and key genes was constructed (Figure 9). Fourteen drugs affected four genes, eight drugs affected three genes, nineteen drugs affected two genes, and thirty-five drugs affected one gene. In particular, among drugs that affected four genes, disulfiram, cannabidiol, and dextroamphetamine have been used to mitigate the cocaine response. Many studies have shown that disulfiram might reduce cocaine use in patients with cocaine dependence (Petrakis et al., 2000; Gaval-Cruz and Weinshenker, 2009; De Mulder and Dom, 2012; Kosten et al., 2013). In the experiment conducted by Petrakis et al. (2000), disulfiram inhibited dopamine β-hydroxylase, resulting in dopamine overdose and decreased norepinephrine synthesis, possibly weakening cocaine cravings, leading to reduced cocaine use. Dextroamphetamine, a central nervous system stimulant, has been found to be a treatment for cocaine dependence (Grabowski et al., 2001; Shearer et al., 2003; Palis et al., 2021; Ndiaye et al., 2022). In experiments on rats by Chiodo and Roberts (2009), sustained dextroamphetamine treatment was found to weaken the potentiating effect of cocaine. Cannabidiol (CBD) is one of the main components of cannabis, and multiple studies have shown that CBD may act as a therapeutic drug for substance abuse (Katsidoni et al., 2013; Calpe-López et al., 2019; Anooshe et al., 2021; Karimi-Haghighi et al., 2022). Recent research showed that CBD can be effective in reducing the reward and reinforcement effects of addictive drugs (Galaj et al., 2020). Among the drugs that affected the three genes, diazepam and melatonin might be useful therapeutic agents for reducing cocaine abuse (Takahashi et al., 2017; Barbosa-Méndez et al., 2021; Sanchez et al., 2022).

In twenty-two drugs affecting the expression of three or four key genes, five (disulfiram, dextroamphetamine, diazepam, cannabidiol, and melatonin) have been validated in the literature to reduce cocaine abuse and be used to treat cocaine addiction. They are distributed among drugs used in addictive diseases, psychoanaleptics, and antiepileptic drugs. We, therefore, speculated that seventeen drugs that affected the expression of three or four genes in these three classes might play the same role.

The Drugbank database was used to analyze the status of these seventeen drugs. All of them have been approved by the FDA for the treatment of other diseases, and their effects on cocaine addiction are still being studied. Twelve of these drugs, namely, disulfiram, nicotine, fluoxetine, donepezil, caffeine, amphetamine, cannabidiol, desipramine, valproic acid, dextroamphetamine, carbamazepine, and methylphenidate, are currently in clinical trials for the treatment of cocaine addiction.

The CTD was then used to analyze the relationship between seventeen drugs and cocaine addiction, and the results showed that disulfiram, dextroamphetamine, caffeine, fluoxetine, methylphenidate, desipramine, scopolamine, valproic acid, diazepam, haloperidol, donepezil, clozapine, carbamazepine, and cannabidiol were chemicals with known or potential therapeutic effects in cocaine addiction. Disulfiram, dextroamphetamine, nicotine, fluoxetine, caffeine, methylphenidate, desipramine, diazepam, scopolamine, amphetamine, haloperidol, and melatonin were chemical substances related to cocaine addiction or may play a role in the etiology of cocaine addiction. Although all seventeen drugs have been confirmed in the CTD to be associated with cocaine addiction, further research is needed to determine whether these drugs can be used to treat cocaine addiction.