Detecting suicide risk in bipolar disorder from lymphoblastoid signatures

Suicide is a deeply complex and tragic phenomenon that continues to be a significant public health concern worldwide. Annually, the global incidence of suicide surpasses 700,000 fatalities exhibiting variations observed among distinct age cohorts and geographical regions [1, 2]. Suicide often unfolds within the framework of psychiatric illness [3]. Research findings show that more than 90% of individuals who succumb to suicide have a diagnosed psychiatric illness and most individuals who attempt suicide also experience such disorders [4, 5]. A significant number of individuals diagnosed with BD attempt suicide at least once in their lifetime, and approximately 20% die by suicide [6, 7]. Overall, BD carries the highest risk of suicidal behavior among all psychiatric disorders [8,9,10]; and currently, it is very difficult to identify those BD patients that are at a high risk of suicide.

The risk associated with suicide in BD varies depending on the disease attributes and the phase or stage of the illness. Existing studies suggest that bipolar disorder type II (BD II) has a higher suicide rate than bipolar disorder type I (BD I). This observation emphasizes the multifaceted severity of BD II [11] with major depressive episodes, carrying the highest risk of suicide among individuals with BD [12]. Following closely are BD patients with mixed episodes [13]. Conversely, during manic episodes, marked by elevated mood and energy levels, the patient has the lowest risk of suicide [14]. Additionally, individuals with rapid cycling BD face an increased risk compared to those without rapid cycling patterns [15]. When considering sex differences certain studies suggest that men exhibit higher rates of fatal suicide attempts, whereas women have more non-fatal suicide attempts [16, 17]. However other research indicates no sex disparity in suicide attempts [18].

In patients with BD, a history of prior suicide attempts emerges as one of the most robust predictors of both fatal and non-fatal suicide attempts [19, 20]. While a significant portion, approximately 56% of those who die by suicide do so on their first attempt [19, 21] the majority of individuals who attempt suicide non-fatally do not ultimately die by suicide [22]. A few studies have highlighted the strong association between a family history of suicide and suicide deaths among individuals with BD [23, 24]. Thus, it confirms the presence of both a family history and a personal history of suicide attempts as a critical suicide risk factor in BD [14, 25]. Other factors include early life trauma such as childhood abuse or stress, which are linked to higher rates of suicide attempts in individuals with BD [26,27,–28]. Moreover, the combination of childhood abuse and drug usage drastically increases the likelihood of suicide attempts [29]. Psychosocial stressors, including interpersonal conflicts and occupational issues, have also been associated with increased suicide risks [30, 31].

Genome-wide association studies (GWAS) have paved the way in identifying genetic loci associated with suicide across psychiatric disorders [32, 33]. Notable genes including 5-HTT, SLC6A4, and the 5-HT1 to 5-HT7 series involved in serotonergic neurotransmission, along with tryptophan hydroxylase genes TPH1 and TPH2, which are identified for their crucial roles in the neurotransmission dysfunctions tied to various suicidal behaviors [34,35,36,37,38,39,40,–41]. Research also extends to genes such as BDNF and its receptor NTRK2, and others like COMT and MAPK1, exploring their links to suicide risk [42,43,44,45,46,–47]. Furthermore, genetic variations in genes like AKT1P, GSK3B, and the Alpha2a-adrenergic receptor (ADRA2) are associated with impulsivity, a significant risk factor in BD suicide [48,49,50,–51]. Additionally, the PENK gene, along with IL-7 and TMX3, involved in stress and anxiety responses, were linked to suicidal behavior in BD patients [52]. Comprehensive multi-ancestry GWAS analysis has pinpointed genes such as DRD2, FURIN, NLGN1, SOX5, PDE4B, and CACNG2 as involved in suicidal behavior, alongside identifying potential biomarkers within the immunoglobulin gene family for BD [53, 54].

Several studies have demonstrated the importance of LCLs as a powerful tool for studying psychiatric disorders, including BD, due to their affordability and scalability [55, 56]. LCLs, derived from patient blood samples, provide a patient-specific model that allows for a large-scale genomic analysis. These advantages have made LCLs a popular choice for investigating the genetic factors associated with BD. For example, in an early study, Morag et al. [57] applied a transcriptome-wide approach using LCLs to identify biomarkers of Selective serotonin reuptake inhibitors (SSRI) response. Their group found significantly lower expression of the CHL1 gene in LCLs with high paroxetine sensitivity, along with differential expression of multiple genes involved in synaptogenesis and psychiatric disorders, demonstrating the utility of LCLs for antidepressant response prediction [57]. Later, Squassina et al. [58] utilized LCLs to identify genetic factors linked to lithium response in BD and found IGF-1 to be significantly over-expressed in lithium responders compared to non-responders [58]. Similarly, Milanesi et al. [59] used LCLs to identify genes like HDGFRP3 and ID2, which showed differential expression between lithium responsive (LR) and non-responsive (NR) BD patients, suggesting their potential as biomarkers for treatment outcomes [59]. In addition, microRNA expression profiling of LCLs from BD patients has revealed specific miRNAs role such as miR-4286, that are differentially expressed in individuals with a history of suicide attempts, highlighting the utility of LCLs for identifying molecular signatures associated with suicide risk [60]. Most recently, Mizrahi et al. [53] from our group identified a set of DEGs in the immunoglobulin gene family from LCLs of BD patients, and trained ML models that successfully predicted lithium response in independent cohorts, demonstrating the scalability and clinical relevance of this approach [53]. These findings have underscored the potential of LCLs in exploring the genetic underpinnings of psychiatric conditions and their potential in identifying biomarkers for personalized treatment strategies.

In this study, we sought biomarkers of suicide risk in BD patients by analyzing transcriptomic differences in LCLs and applying ML algorithms to classify high vs. low risk based on gene expression profiles. To conduct this study, we extracted RNA from LCLs obtained from patients who had been followed and monitored prospectively for several years. Among these patients, some eventually died by suicide, while others with BD never attempted suicide. We also retained a third group of patients to assess the final performance of our classifier after conducting train-test cross validations using the first two groups. Interestingly, when analyzing biological pathway enrichment between the ‘SUICIDE’ and ‘NON-SUICIDE’ groups, we identified pathways related to brain and psychiatric diseases, despite the RNA source being LCLs. The patients were categorized into two groups: those who died by suicide, classified as high risk, and those who did not, considered low risk, based on long-term monitoring.

Ethics statement and LCLs collection

The study received approval from the University of Cagliari’s local Ethics Committee (approval number: 348/FC/2013), and the Research Ethics Committee at the Nova Scotia Health Authority approved the study in Halifax. All participants provided informed written consent at the time of recruitment, following a thorough explanation of the study protocol. The procedures adhered to the ethical principles outlined in the Helsinki Declaration of 1975.

Data set

Diagnoses for all participants were established according to the Research Diagnostic Criteria [61] and DSM-IV guidelines. These diagnoses were based on semi-structured personal interviews (Schedule for Affective Disorder and Schizophrenia Lifetime Version (SADS-L)) [62, 63] and thorough reviews of the patients’ medical histories. Interviews were also conducted with first- and second-degree relatives, using SADS-L and diagnoses were assigned in accordance with RDC and DSM-IV standards. For relatives who could not be interviewed in person, diagnoses were made using the Family History-Research Diagnostic Criteria, which involved gathering information from at least two informants and reviewing available medical records.

In this study, we conducted RNA sequencing on samples from 20 patients diagnosed with BD. All patients included in the study are of Caucasian ethnicity. They were treated with various medications, including Benzodiazepines (BZDs), Carbolithium (CAR), Carbamazepine (CBZ), Clonazepam (CLZ), Promazine (PMZ), Flurazepam (FLZ), Alprazolam (ALP), Quetiapine (QTP), Lansoprazole (LZP), Lamotrigine (LTG), Venlafaxine (VFX), Bupropion (BPR), Risperidone (RSP), Fluoxetine (FLX), Olanzapine (OLZ), Etizolam (ETZ), Haloperidol (HAL), Amitriptyline (AMT), Lorazepam (LRZ) and Levothyroxine (LT4). The specific treatment regimen for each individual is detailed in Table 1. In this study, all 20 patients were divided into two main groups. The first 13 patients were used for training and validation. This subset included a “SUICIDE” group (n = 6), comprising patients who died by suicide, and a “NON-SUICIDE” group (n = 7), consisting of individuals who were prospectively monitored for a median of 4 years (range: 1–7 years) without any suicide attempts or a known family history of suicide. These 13 patients were randomly divided into training and validation sets. The ML algorithm was evaluated on the validation set (unseen by the algorithm), which included both high-risk (died by suicide) and low-risk patients. Detailed information on all 20 patients is provided in Table 1. The remaining 7 patients (apart from the 13 used for training and validation) were used exclusively for testing. These patients were not included in the training or validation phases and were therefore entirely unseen by the algorithm. This group had mixed risk profiles (e.g., family history or non-fatal suicide attempts) and was used to preliminarily assess the algorithm’s ability to generalize and classify new, unobserved data. Additional information about the samples will be available in the GitHub file at this link (https://github.com/omveersharmanet/Predicting-Suicide-Risk-in-Bipolar-Disorder patients.git).

LCL growth and expansion

We received LCLs derived from Peripheral blood mononuclear cells (PBMCs) from both groups of patients. Briefly, PBMCs were isolated from blood samples using BD vacutainer according to the manufacturer’s guidelines. Following which they were exposed to Epstein-Barr virus (EBV) harvested from the B95-8 cell line. Within a week post-infection, the cells transformed and began to aggregate, forming large clusters over time. The LCLs were cultured in T25 tissue culture flasks, using a complete RPMI medium composed of RPMI 1640 (Biological Industries, Cat no: 01-100-1A), 1X Anti-Anti (Thermofisher Scientific, Cat no: 15240062), 1% Glutamax (Thermofisher Scientific, Cat no: 35050061), 1% Sodium pyruvate (Thermofisher Scientific, 11360070), and 15% heat-inactivated fetal bovine serum (Sigma, Cat no: F9665), as previously outlined [64]. Media changes were performed every other day. Passaging of the LCL was performed to maintain optimal cell growth and prevent over-confluency. Cells were sub-cultured when reaching a density of 200,000 cells/ml to sustain a logarithmic growth.

RNA extraction, sequencing, and analyses

Total cellular RNA was extracted from 5 million LCLs using 800 μl of TRIzolTM reagent (Thermofisher Scientific, Cat no. 15596026) and subsequently chilled on ice for 5 min before long-term storage at −80° Celsius. The total RNA was carefully isolated using the zymo RNA clean& concentrator kit, according to manufacturer’s protocol. Following the extraction, the quality and integrity of the RNA samples were quantified using ND-1000 Nanodrop spectrophotometer (Thermo Sci). All RNA samples displayed an RNA integrity number (RIN) within 7–8, confirming both good quality and suitability for sequencing.

RNA sequencing and analysis Libraries of RNAs extracted from LCLs derived from BD patients were prepared using the TruSeq RNA Library Prep Kit v2 (Illumina), adhering to the manufacturer’s guidelines. Quality control was conducted on the raw FASTQ files using FastQC (v0.11.5) [65]. Reads were aligned to the human genome (GRCh38.104) and quantified with STAR (v2.7.9a) [66]. Differential gene expression analysis was carried out using DESeq2 (v1.34.0) [67]. To mitigate false positives in identifying DEGs, a false discovery rate (FDR) analysis was applied using the Benjamini-Hochberg (BH) procedure to adjust p-values for multiple hypothesis testing [68]. Genes were considered significant DEGs if they met the criteria of a p-value < 0.05 and a log2 fold-change ≥ |0.58|, with a total of 841 genes meeting these thresholds. The step-by-step processStructural overview of the integrated RNA-seq data analysis and machine learning framework. Structural overview of the integrated RNA-seq data analysis and machine learning framework. for identifying DEGs is illustrated in Fig. 1a. The top 50 genes were selected as predictors for further machine-learning analysis.

Fig. 1: Structural overview of the integrated RNA-seq data analysis and machine learning framework.

a A Flowchart illustrating the process of identifying differentially expressed genes (DEGs) and training and testing the ML models. b A step-by-step training and testing of the ML models.

Machine learning predictor analysis

We assessed five distinct supervised classifiers for their ability to predict suicidal and non-suicidal outcomes; Logistic Regression (LOR) is a fundamental statistical model utilized for binary classification, leveraging a logistic function to estimate probabilities and optimize parameters through maximizing data likelihood, though its performance may falter with non-linear data relationships. The K-Nearest Neighbors (K-NN) algorithm, a non-parametric method, classifies new instances by using the majority vote from the ’K’ nearest training samples, but struggles with large or high-dimensional datasets due to the heavy computational demand of distance calculations. Support Vector Machine (SVM) creates a hyperplane in high-dimensional space to separate classes with optimal margins. It adapts well to nonlinear boundaries using the kernel trick, although tuning its hyperparameters can be challenging and resource-intensive [69, 70]. Random Forest (RF) enhances robustness and accuracy by building multiple decision trees and integrating their outcomes [71], which are suitable for various tasks but potentially complex and difficult to interpret. In our previous research, we discovered that SVM and RF were the most effective machine-learning algorithms in distinguishing between lithium responders and non-responders among BD patients [72,73,74,–75]. Inspired by the human brain, neural Networks (NN) excel in pattern recognition and learning complex nonlinear relationships, making them powerful yet resource-heavy and challenging to manage due to their propensity for over-fitting [76]. Each algorithm offers unique strengths and limitations, with the choice often depending on data characteristics and specific problem needs. Each of these classifiers has been carefully selected to address different aspects of the predictor analysis in our study, taking into account their respective strengths and limitations according to the specific requirements of the data and the task at hand. The LOR, RF, and NN classifiers gave the best performance with high accuracy (>96%). As mentioned earlier, the ML algorithms were evaluated on the validation set (unseen by the algorithm), which included both high-risk (died by suicide) and low-risk patients. The models achieved high accuracy on this validation set, supporting its ability to distinguish between risk categories based on known outcomes.

The feature selection process (predictors/genes) is illustrated in Fig. 1b. We have employed an intensive approach that combines leave-one-out cross-validation (LOOCV) with an exhaustive feature selection process to identify the most crucial features for our model, especially useful for smaller dataset. Specifically, we utilized LOOCV for our dataset, which consists of 13 samples. This method involved training our model on 12 samples and testing on the remaining one, cycling through iteratively so that each sample served as the test set once, resulting in 13 total iterations. Concurrently, during each LOOCV iteration, we examined all possible combinations of up to 10 genes from the initial top 50 DEGs, assessing the LR model’s accuracy with these subsets on the left-out sample, where the outcome was binary—100% or 0% accuracy based on correct or incorrect predictions. Successful combinations led to a tally of each gene’s appearance frequency, and after 13 iterations, the top 20 genes were identified based on this frequency. This process was repeated to yield 13 lists of the top 20 genes, from which we further distilled the data to pinpoint the top 10 consistently appearing genes. Finally, to evaluate our ML models, these top 10 genes were used to split the full dataset into training and testing sets multiple times (50 splits with a 50–50% ratio), and the mean accuracy from these tests provided a robust assessment of the model’s effectiveness. Kindly note that the top 30 and top 50 genes were selected based on the lowest adjusted p-values from differential expression analysis, reflecting strong statistical significance. In contrast, the top 10 genes identified by the ML algorithm were chosen from the top 50 to maximize feature diversity. As p-values capture statistical differences and ML emphasizes predictive patterns, overlap between the two gene sets is not expected. The DEGs, top 50, top 30, and ML-selected top 10 genes are provided in Supplementary Table 1.

The effectiveness of the approach was validated by testing five supervised classification algorithms using the specified input features for binary classification. The outcomes of the model predictions were classified into four categories: true positives (t_p), false positives (f_p), true negatives (t_n), and false negatives (f_n). Evaluation metrics included accuracy, the Receiver Operating Characteristic (ROC) curve, and the confusion matrix. To ensure robustness, cross-validation methods were applied, with results being averaged or aggregated for the confusion matrix across iterations. Performance was assessed using the Area Under the Curve (AUC) of the ROC and metrics from the confusion matrix such as Accuracy and Precision. Figure 1b illustrates the feature selection process (predictors/genes) and the training and testing of ML models using the selected genes. Accuracy of ROC:

Where, y′ represents the model’s prediction, while y denotes the actual class.

Our research focused on developing biomarkers that can be relatively easy to extract from biological samples (as compared to neurophysiology from patients’ neurons) [53, 72, 77]. For this, we have grown LCLs of BD patients who died by suicide (6 samples) and BD patients who have been followed over many years and have not attempted suicide, nor have any family history of suicide. For this, we performed the following steps:

1. (1)

   Identifying DEGs between patients who died by suicide (‘SUICIDE’) compared to those that are at a low risk (no attempts over many years and no family history ‘NON-SUICIDE’),

2. (2)

   Selecting the most predictive genes, and

3. (3)

   Training and testing ML models first on these two groups and then on another set of patients who are at an increased risk of suicide, as shown in Fig. 1b.

The identification of DEGs in LCLs of BD patients (‘SUICIDE’ and ‘NON-SUICIDE’ Groups)

To identify biomarkers for suicide risk, we performed sequencing of RNA extracted from LCLs of 20 BD patients. Of these, six patients died by suicide (‘SUICIDE’), seven patients have never attempted suicide and do not have a family history of suicide (‘NON-SUICIDE’), and the remaining seven patients, who have a family history of suicide and some have attempted non-fatal suicide, are classified as ‘test (unseen samples)’. All patients were followed over the years. Differential gene expression analysis was conducted between the first two groups of BD subjects (‘SUICIDE’ and ‘NON-SUICIDE’), as detailed in the Methods section [2] and in Fig. 2a. We identified 841 DEGs with a fold change greater than 1.5 (log2 fold-change ≥ |0.58|) and a p-value of less than 0.05. These are presented in Fig. 2b as a heat-map. Figure 2b displays the gene counts for the 30 most significantly DEGs based on the lowest p-values. DEGs between ‘SUICIDE’ and ‘NON-SUICIDE’ LCLs exhibited enrichment in pathways related to “Primary immunodeficiency” such as LCK, CD8A, AICDA, ZAP70, and CIITA. Specifically, LCK, AICDA, ZAP70, and CIITA demonstrate higher gene counts in ‘SUICIDE’ LCLs compared to ‘NONSUICIDE’ LCLs, indicating increased expression levels associated with suicidal behaviour. However, CD8A shows a contrasting pattern, with higher gene counts observed in ‘NON-SUICIDE’ LCLs compared to ‘SUICIDE’ LCLs, suggesting a distinct expression profile in non-suicidal individuals. DEGs between the ‘SUICIDE’ and ‘NON-SUICIDE’ groups were enriched in Immunoglobulin genes, such as IGHV1-3, IGHV4-31, IGHV3-33, IGHV1-46, IGHV4-59, IGHV4-61 were enriched in ‘NON-SUICIDE’ samples whereas IGKV4-1 was enriched in ‘SUICIDE’ samples. 0.8% (470 genes) of the total human genes are immunoglobulin genes. However, among the DEGs, 2.3% (19 out of 841 DEGs) are immunoglobulin genes. When examining DEGs with minimum counts of 10, 100, and 1000, the proportion of immunoglobulin-related genes increases to 2.3, 3.1 and 5.2%, respectively. This pattern indicates a significant enrichment of immunoglobulin-related genes in the DEGs, as illustrated in Supplementary Fig. 1(a). Furthermore, we also compared our DEGs with established GWAS datasets along with recently published review articles [39,40,–41, 78]. A broader view of the GWAS dataset and its overlap with our gene sets is illustrated by venn diagram in Supplementary Figs. 2 and 3.

Fig. 2: Flow chart, Genes and pathways that distinguish the two groups of ‘SUICIDE’ and ‘NON-SUICIDE’.

a Flow chart (b) A heatmap of log10 gene count of the DEGS between the two groups of ‘SUICIDE’ and ‘NON-SUICIDE’. The heat map was organized based on the descending order of the subtraction of the average read count of the two groups (‘SUICIDE’-‘NON-SUICIDE’). c The gene count of the top 30 genes with the lowest p-values is displayed, with red indicating suicidal individuals and blue indicating non-suicidal individuals.

Associated dysregulated pathways and diseases in LCLs of BD patients (‘SUICIDE’ vs. ‘NON-SUICIDE’)

Further analysis of the DEGs included exploring KEGG pathways, disease associations, and functional annotations, as depicted in Fig. 3. The KEGG pathway analysis highlights the intricate biochemical and physiological networks potentially involved in the underlying conditions studied, as shown in Fig. 3a. For example, pathways such as primary immunodeficiency, phenylalanine, tyrosine, tryptophan biosynthesis, and hematopoietic lineage are significantly dysregulated between BD patients from the two groups (‘SUICIDE’ vs. ‘NON-SUICIDE’). Additionally, the Intestinal Immune Network for IgA Production and Primary Immunodeficiency pathways suggest a significant role of immune system functions in the pathology, potentially linking systemic immune responses to neurological and psychiatric conditions.

Fig. 3: Dysregulated pathways, disease, and functional annotations in LCLs of BD patients who died by suicide compared patients at a low risk of suicide.

a A diagram displaying the dysregulated KEGG pathways associated with the DEGs between the two groups. b The most prevalent conditions in the DisGeNET database, such as Libman-Sacks Disease, Systemic Lupus Erythematosus, Thrombosis, Neoplasm Invasiveness, Gliosis, Astrocytosis, and various conditions related to pregnancy loss and disorders. c Prominent GAD disease classes, including Cardiovascular, Hematological, Immune, Psychiatric, Renal, and others. d Significant biological processes like Cell Adhesion, Immunity, Chemotaxis, and Inflammatory Response. e Molecular functions associated with the DEGs, including Ion Channels and Growth Factors.

In the Disease DisGeNET analysis, a wide array of conditions linked to the DEGs was identified, indicating a broad impact of these genetic variations across various diseases, as shown Fig. 3b. The list includes autoimmune and inflammatory diseases such as Libman-Sacks Disease and Systemic Lupus Erythematosus, as well as conditions related to abnormal blood clotting, like thrombosis. Importantly and surprisingly, many brain-related pathways came up in the analysis, despite the origin of the cells (LCLs). These pathways include BD Schizophrenia, Cerebral infraction Left and Right hemispheres, Subcortical infraction, Autistic Disorder, neuroinflammatory diseases like Gliosis, Astrocytosis, and more.

The Disease GAD (Genetic Association Database) list encompasses a wide range of disorders, highlighting the complex interplay between genetics and various health conditions, as shown in Supplementary Fig. 1(b). Diseases such as Hypertension, Asthma, and Coronary Artery Disease reflect chronic conditions that significantly impact public health and are influenced by both genetic predispositions and lifestyle factors. Autoimmune and inflammatory diseases like Systemic Lupus Erythematosus and Rheumatoid Arthritis, as well as metabolic disorders such as Type-2 Diabetes, underscore the genetic basis of immune system dysfunction and metabolic dysregulation. Furthermore, the inclusion of conditions like schizophrenia and substance use disorder illustrates the genetic factors contributing to psychiatric disorders and addictive behaviors, demonstrating the broad spectrum of diseases that genetic research can potentially address. In our research, we conducted an in-depth analysis of diseases characterized by specific up-regulated genetic or molecular markers, focusing on how these enhancements influence disease pathology, as shown in Supplementary Fig. 1(b). This exploration has been pivotal in identifying key biological pathways and potential therapeutic targets, enhancing our understanding of their underlying mechanisms.

The GAD Disease class encompasses a broad spectrum of medical categories, reflecting the diverse genetic underpinnings and environmental influences on various health conditions, as shown in Fig. 3c. Interestingly, our data indicates that neurological and developmental pathways are dysregulated even in the white blood cells of the patients. Another interesting pathway that was dysregulated was “AGING”, which aligns with previous reports of accelerated aging in BD patients having suicide attempts/ideation [79,80,–81]. “PSYCH” was one of the top dysregulated pathways in the LCLs. The “IMMUNE” pathway is among the top 3 dysregulated pathways. Immune pathways have been shown repeatedly to be dysregulated in BD [82]. The most dysregulated pathway was “CARDIOVASCULAR”. Associations between BD and heart defects were previously shown [83]. Our data shows that these pathways are even more dysregulated in BD patients, who are at a high risk of suicide.

We also conducted a detailed analysis across five functional annotation categories for which the DEGs between the groups were enriched including UP KW BIOLOGICAL PROCESS, UP KW CELLULAR COMPONENT, UP KW MOLECULAR FUNCTION, UP KW PTM, and UP SEQ FEATURE as shown in Fig. 3d, e and Supplementary Fig. 1(c)–(e). Specifically, within the UP KW BIOLOGICAL PROCESS category, we identified significant dysregulation in processes such as cell adhesion, which is vital for cellular assembly and integrity, and various immune-related functions, including adaptive and innate immunity, underscoring the genes’ pivotal roles in systemic defense mechanisms, as shown in Fig. 3d. Additionally, this category highlighted the involvement of genes in essential developmental and reparative processes such as angiogenesis and osteogenesis, crucial for maintaining vascular and skeletal health. In our study, the UP KW MOLECULAR FUNCTION revealed critical roles of various proteins in cellular mechanisms, with a particular emphasis on ion channels, sodium channels, and voltage-gated channels, which are essential for electrical signaling and cellular communication and have been shown to alter in neurons derived from BD patients [84], as shown in Fig. 3e. Proteins such as cytokines and growth factors were also identified as changed in LCLs of suicide victims [85, 86]. These proteins are important for cell signaling and regulatory processes that influence cell growth, differentiation, and immune responses, suggesting implications for suicide [87, 88]. Additionally, the identification of host cell receptors for virus entry, integrins, and receptors underscores the relevance of cell-environment interactions in suicide, which facilitates crucial processes such as viral entry, cellular adhesion, and signal transduction across cellular membranes.

Selecting the best predictor genes

We initiated this study with the goal of developing biomarkers that could be quickly and easily implemented in clinical settings. Consequently, it was essential to establish a robust protocol that facilitates the prediction of BD patients who are at a high risk of suicide using RNA-seq datasets, as illustrated in Fig. 2a. As described in the methodology section [2], the top 50 genes were selected from a pool of 841 DEGs to serve as predictors for ML algorithms. We implemented a rigorous approach to identify the most dominant and predictive genes using LOOCV and exhaustive feature selection, well-suited for our dataset. Our method involved training the model on 12 samples and testing on the remaining one, cycling through each sample as the test set, completing 13 iterations. During each iteration, we tested all combinations of up to 10 genes from the initial 50, assessing the LOR model’s accuracy on the left-out sample, which was binary either 100 or 0% based on prediction accuracy. This allowed us to identify the top 20 genes appearing most frequently across iterations. After repeating this process, we narrowed down to the top 50 genes consistently identified across the 13 lists. The gene counts and heat maps for these top 50 genes are shown in Fig. 4a, b. We further narrowed down the top 10 genes to train and test all ML models, as shown in Fig. 1b.

Fig. 4: Top predictive genes for suicide risk.

a A heat map, the heat map was organized based on the descending order of the subtraction of the average read count of the two groups (‘SUICIDE’-‘NON-SUICIDE’). b Gene counts for the most dominant predictor for ML.

Training and testing of machine-learning models

After identifying the top 10 genes (LYL1, LMO7, DDAH2, PAG1, LCK, AP001372.2, INPP5D, SNX8, TBC1D4, GJC1), we proceeded to train and test various ML models. We employed five types of supervised classification algorithms: LOR, RF, K-NN, SVM, and NN. To avoid over-fitting, we used cross-validation methods, as detailed in the Methods section [2] (refer to Fig. 1b). We trained and tested the models using subsets varying from 1–10 predictors, as illustrated in Fig. 5a. Notably, we found that 8 features provided high accuracy of more than 95%.

Fig. 5: Prediction accuracy by varying the number of genes.

a The variation in model accuracy for a LOR model in relation to the number of genes used as predictor features, ranging from 1–10 genes. b The accuracy of various ML models when specifically using 8 predictors.

We narrowed down to 8 genes with high predictive power of suicidal and non-suicidal outcomes, as indicated in Fig. 5b. Here, we divided the data with 50% for training and 50% for testing, a process repeated 50 times to ensure the robustness of our ML algorithms. The mean classification accuracies with standard deviations across repetitions were: LOR: 99.67% ± 2.36%, RF: 95.67% ± 12.05%, K-NN: 66.33% ± 29.45%, SVM: 88.00% ± 19.06%, and NN: 99.00% ± 4.00%. Both LOR and NN showed higher accuracy rates compared to the other ML models. Figure 6 displays the performance of these classifiers through the area under the ROC curves and confusion matrices, with the following AUC scores: LOR: 1.00 ± 0.00, Random Forest: 0.99 ± 0.00, K-NN: 0.69 ± 0.08, SVM: 0.91 ± 0.05, and NN: 1.00 ± 0.00. The ML algorithms demonstrated high accuracy on the validation set, confirming its effectiveness in distinguishing between risk categories based on known outcomes.

Fig. 6: ROC curves and confusion matrices for machine learning models.

a LOR, (b) SVM, (c) RF, (d) NN, and (e) K-NN, along with (f) an overlay of all ROC curves (NN and LOR exhibit overlapping ROC curves).

Furthermore, the remaining 7 patients (excluding the 13 used for training and validation from the total of 20 patients) were reserved exclusively for testing. These individuals were not involved in any phase of training or validation and were therefore entirely unseen by the algorithms. They have mixed risk profiles (e.g., family history or non-fatal suicide attempts) and were used to assess the algorithm’s ability to generalize and accurately classify previously unobserved data. However, it is important to note that psychosocial factors such as sex, age, BD subtype, symptomatic stage and psychiatric comorbidities, or family history were not included as input features for the ML algorithms. On the other hand, family history and the patient’s suicide history were analyzed separately to examine its correlation with the algorithm’s predicted suicide risk.

The trained models were also used to assess the status of 7 BD patients whose conditions were not previously identified, yet some of their clinical data suggests that they are at a high risk of suicide, such as a first-degree relative who died of suicide and previous suicide attempts. For this classification, we took the majority decision by each of the classifiers (‘Combined Results’). The analysis revealed that patients UOHL027 and UOHL031 are at a low risk of suicide, whereas the other five BD patients are classified as at a high risk. These findings are detailed in Table 2. Interestingly, the BD patient with no previous suicide attempts and no family history was classified as a low risk for suicide. The other one that was classified as a low risk had no previous attempts and one family member who died by suicide. The rest of the patients were classified as a high-risk for suicide. These patients had previously attempted suicide or had a few family members who died by suicide (with one patient UOHL030 having just one family member who died by suicide). The test sample was not ideal, since it did not comprise of patients who dies by suicide, but rather high-risk or low-risk patients based on their family and suicide attempts history. This may be a limitation of the current approach, however, live cell samples of patients who died suicide are naturally extremely hard to obtain. Because of the rarity and complexity of suicide events, and their inherent unpredictability, obtaining large, balanced datasets with well-defined outcomes remains a significant challenge in this field. This cohort was acquired after long years of following BD patients at the clinic.

In this study, we identified distinct gene expression signatures associated with increased suicide risk by utilizing RNA sequencing data from LCLs of BD patients who died by suicide, compared to those at low risk of suicide. Notably, the genes involved in cardiovascular, ion channel, and cell adhesion pathways. After identifying genes and pathways linked with suicide, we have further used these as features to train classification algorithms in order to construct a predictor for suicide risk in BD using biomarkers that are easily and cheaply available without much stress to the patient. We initially identified 841 genes that were differentially expressed in LCLs of BD patients who died by suicide compared to patients at a low risk of suicide. The low-risk patients were followed and monitored over the years and have not attempted suicide, nor did they have a family record of suicide. We used several pathways, disease, and functional annotation analyses to first identify possible mechanisms of suicide. Importantly, our biological samples were LCLs immune-related cells. We therefore expected to find immune-related changes and did not expect to find brain-related implications. For brain-related mechanisms, we assumed that induced pluripotent stem cell (iPSC) work related to differentiating the cells into neurons would be required. We were very surprised, therefore, to find many neuronal pathways that are altered between patients at a high risk of suicide and patients with the low risk. This suggests that when differentiating neurons, we will find neurophysiological changes in neurons derived from patients at a high risk of suicide compared to those at a low risk. It is worth noting that numerous studies have employed diverse cell line models in investigating BD and suicide risk, reflecting the complex nature of these conditions. A detailed summary on recent findings related to BD and suicide from 2020–2025 are provided in Supplementary Table 2 [60, 81, 89,90,91,92,93,94,95,96,97,98,99].

Among the primary findings, DEGs were enriched for the primary immunodeficiency pathway between LCL from BD patients who died by suicide and patients who are at a low risk of suicide. For example, genes such as LCK (log2 fold change of 0.81), AICDA, (log2 fold change of 1.329), and CIITA (log2 fold change of 0.975) exhibited increased expression levels in LCL from BD patients who died by suicide compared to those at a low risk of suicide. These genes are critical for immune system functioning; LCK is involved in T-cell receptor signaling [100]. CIITA (Class II Major Histocompatibility Complex Transactivator) regulates MHC (Major Histocompatibility complex) class II genes [101]. AICDA is crucial for B-cell antibody response [102]. The differential expression of these genes might reflect an underlying dysregulation in the immune response, which could be contributing to the pathophysiology of suicide. On the other hand, the CD8A gene (log2 fold change of 2.008), which is vital for T-cell mediated cytotoxicity [103], showed lower expression in LCLs from people who died by suicide, suggesting a potential protective role of this gene against suicide. This pattern of gene expression adds an intriguing layer to our understanding of how genes associated with both high and low risk can be linked to immune dysregulation and psychiatric conditions, particularly in the context of suicide in BD. Continued surveillance of these genes is crucial, as it will enable further exploration into the complex relationships between immune function disturbances and suicide.

Further supporting the link between immune function and the brain, our study also highlighted significant alterations in cell adhesion and ion channel pathways in LCLs of BD patients who died by suicide compared to the patients at a low risk. These pathways are known to influence neuronal connectivity and excitability [104, 105]. From the cell adhesion pathway, genes like NTNG1 and NRXN1 showed increased expression, with log2 fold changes of 2.27 and 1.42, respectively in LCLs of BD patients who died by suicide compared to BD patients with a low risk of suicide. Both genes are crucial for axonal guidance and synaptic stability [106, 107]. Additionally, CD8A was upregulated with a log2 fold change of 2.008 and. This gene also participates in the primary immunodeficiency pathway. The gene CLDN10 is a key component of the blood-brain barrier (BBB) tight junctions [108] and was down-regulated with a log2 fold change of −2.16 in LCLs of people who died by suicide. Dysregulation in cell adhesion genes can lead to altered excitatory/inhibitory balance within the neural circuit [109,110,111,–112].

Ion channels are crucial in modulating neuronal excitability and signal transduction, with growing evidence linking their dysregulation to psychiatric disorders [113, 114]. Our study highlights significant gene expression differences in BD patients who died by suicide compared to those with a low risk of suicide. Specifically, voltage-dependent calcium channel genes such as CACNB2 (Calcium Voltage-Gated Channel Auxiliary Sub- Unit Beta 2) showed a notable increase (log2 fold change of 1.318) in BD patients who died by suicide. In contrast, CACNA1A (Calcium Voltage-Gated Channel Sub-unit Alpha1 A) exhibited a log2 fold change of 0.955, predominantly in patients with a low risk of suicide. These genes mediate calcium influx into a wide variety of electrically excitable cells, including cardiac, muscle, neurons, and sensory cells [115]. In a similar pattern, SCN11A and SCN1B genes, which encode for voltage-gated sodium channels [116], showed differential expression. SCN11A (sodium voltage-gated channel alpha subunit 11) was up-regulated (log2 fold change of 1.355) in BD patients who died by suicide, while SCN1B (sodium voltage-gated channel beta subunit (1) showed a higher expression (log2 fold change of 2.33) in patients with a low risk of suicide. This differential expression of calcium and sodium channel genes suggests distinct roles, potentially categorizing them as either low or high predictors of suicide risk. Moreover, genes encoding for calcium-activated potassium channels such as KCNN2 (Potassium Calcium-Activated Channel Subfamily N Member (2) and KCNMB1 (Potassium Calcium-Activated Channel Subfamily M Regulatory Beta Subunit 1) were also highly expressed in BD patients who died by suicide [117] with a log2 fold change of 3.228 and 1.988 respectively. Additionally, the serotonin receptor gene HTR3B also showed a significant up-regulation (log2 fold change of 3.4) in BD patients who died by suicide, reaffirming the well-established link between serotonin signaling disruptions and both depressive and suicidal behaviors [34, 118,119,120]. These insights reveal a complex interplay between genetic factors and ion channel function, further elucidating the multifaceted nature of suicide risk in BD patients.

Likewise, our findings support the established link between ion channel dysfunctions and their role in both cardiovascular and neurological disorders, highlighting their critical roles in cellular excitability and signaling [121]. This connection has been well-documented, with ion channel abnormalities known to contribute to conditions such as cardiac arrhythmias and various neurological maladies, suggesting a shared pathophysiological foundation [122]. In our study, we observed that genes such as KCNN2 (log2 fold change of 3.228), KCNMB1 (log2 fold change of 1.98), SCN11A (log2 fold change of 1.35), CACNB2 (log2 fold change of 1.318) and glutamate ionotropic receptor AMPA type subunit 1 also known as GRIA1 (log2 fold change of 2.85), showed increased expression in LCLs from BD patients who died by suicide, while a few ion channels like SCN1B (log2 fold change of 2.33) and KCNE4 (Potassium Voltage-gated channel subfamily E member 4) (log2 fold change of 4.06) were more prominently expressed in BD patients with a low risk of suicide. This pattern of differential gene expression underscores the significant role these channels play in influencing both neuronal and cardiovascular functions, which in turn affect mood and behavior. The data points to a complex interaction between genetic factors that impact both neurological health and cardiovascular stability, reinforcing the potential of these genes as dual biomarkers for assessing the risk of suicide as well as cardiovascular diseases.

We also identified genes associated with astrocytosis and gliosis. Astrocytosis and gliosis involve the proliferation and hypertrophy of astrocytes, usually in response to central nervous system pathologies like neurodegenerative disease, tumor growth, and trauma [123]. We found five genes to be differentially expressed, among which BDNF and ITGA5 were notably upregulated in BD patients who are at a low risk of suicide, with log2 fold changes of 2.04 and 1.30, respectively. BDNF is pivotal for neuronal survival and development, influencing the growth, differentiation, and synaptic plasticity of neurons [124]. It also plays a critical role in neuronal and synaptic regeneration and repair following injury [125]. On the other hand, ITGA5 is essential for cell adhesion to the extracellular matrix, supporting cell migration during development and wound healing, and is involved in angiogenesis [126, 127]. The roles of neurotrophins like BDNF and integrins like ITGA5 are crucial in guiding neuron positioning, growth, and maturation, highlighting why their low expression might be observed in suicide samples [128, 129]. However, the GRM8 gene was found to be elevated in LCLS of BD patients who died by suicide (log2 fold change of 1.38) and encodes for glutamate metabotropic receptor 8 and plays a significant role in glutamate neurotransmission. This gene’s function is particularly pertinent in BD, where glutamate system dysregulation is often implicated [130]. Additionally, the GRM8 gene’s involvement in processes like astrocytosis can be linked to altered glutamatergic signaling which may influence the proliferation and hypertrophy of astrocytes during neuroinflammatory responses.

Intriguingly, beyond their hematopoietic lineage, LCLs exhibit a remarkable expression profile reminiscent of neuronal activity, as evidenced by RNA sequencing studies [131, 132]. This unexpected convergence underscores the potential of LCLs as a valuable model for probing the genetic underpinnings across the spectrum of brain disorders [133, 134]. In our analysis, several genes, including SHANK2 (log2 fold change 1.96), CACNG5 (log2 fold change 4.93), ANK3 (log2 fold change 2.04), and NRG3, (log2 fold change 4.29) were found to be significantly upregulated in LCLs of BD patients who died by suicide compared to patients who are at a low risk of suicide. SHANK2 plays a crucial role in synaptic transmission and has been associated with autism spectrum disorders and schizophrenia [135, 136]. CACNG5 encodes a subunit of the voltage-dependent calcium channels [137, 138] and has been previously known to be linked to schizophrenia [138]. ANK3 is crucial for encoding ankyrin-G, a protein involved in the structural stability of neurons by anchoring integral membrane proteins to the cytoskeleton [139]. Previous studies have shown that variations in this gene are associated with both BD and schizophrenia [140]. NRG3, encoding neuregulin 3, is a part of the neuregulin family. Neuregulins act through binding to the ERBB family initiating important signaling pathways that influence cell growth, and differentiation [141]. Genetic studies have linked variations in the NRG3 gene to schizophrenia and BD [142]. This concurrent association underscores the potential role of these genes as high-risk factors for suicide in BD. The upregulation of these genes not only reinforces their critical function in neuropsychiatric pathologies but also highlights them as promising biomarkers for assessing suicide risk.

Interestingly, within our set of DEGs, we identified several genes that have previously been associated with lithium responsiveness (LR) and non-responsiveness (NR) in BD. Building upon earlier studies that have already characterized LR and NR genes [75, 143], we performed a comparative analysis to find an overlap between the current DEGs and previously reported gene sets. Notably, we identified 10 common genes between our DEGs and the 229 DEGs previously reported as differentially expressed between LR and control samples [75]. Similarly, we found 120 common genes between our DEGs and the 2764 DEGs reported between NR and control samples, as shown in Supplementary Fig. 4. Additionally, an overlap analysis between our DEGs and the combined set of 2876 DEGs from LR vs. Control and NR vs. Control comparisons revealed 123 genes that were commonly detected in our dataset. Notably, these findings reveal that very few LR genes were detected in our analysis. This observation aligns remarkably well with the clinical characteristics of our study cohort. As shown in Table 1, the initial 13 patients listed in our clinical and demographic data were either classified as NR or did not belong to the LR group. This convergence between our analysis and clinical phenotypes suggests that the predominance of NR gene expression may reflect the underlying biological mechanisms in patients who do not respond favorably to lithium treatment. These findings provide important new evidence that connects LR/NR (DEGs) genes with vulnerability to BD and suicide risk emphasizing the need for further research into these molecular pathways to better understand and mitigate suicide risk in this population. Furthermore, we acknowledge that sex differences, particularly the higher prevalence of suicide among males, may confound gene expression analyses, as certain biological pathways such as those related to cardiovascular function can vary by sex. To investigate this potential influence, we performed principal component analysis (PCA) on the top 20 differentially expressed genes, as shown in Supplementary Fig. 5. The PCA results showed a clear separation between the suicide and non-suicide groups, indicating distinct gene expression profiles. Interestingly, two male patients in the non-suicide group (UOHL020 and UOHL021) clustered more closely with the suicide group, suggesting that these individuals may exhibit transcriptomic signatures associated with increased suicide risk. In contrast, one female patient in the suicide group clustered nearer to the non-suicide group. These findings suggest the possibility of sex-related differences in gene expression associated with suicide risk. However, due to the small sample size, we caution against overinterpreting these trends.

To delve deeper, we enhanced our understanding of the genetic factors that influence BD, particularly in distinguishing between ‘SUICIDE’ and ‘NON-SUICIDE’ tendencies, by refining our research methodology through a robust protocol. Out of 841 DEGs, we selected 50 potential biomarkers to be analyzed via ML techniques. We established a list of the top 10 markers, which showcased consistently high predictive accuracy. These genes are LYL1, LMO7, DDAH2, PAG1, LCK, AP001372.2, INPP5D, SNX8, TBC1D4, and GJC1. Importantly, these top 10 genes differ from the initial 841 DEGs, showing clearer gene count distinctions between ‘SUICIDE’ and ‘NON-SUICIDE’ cases. Our methodology involved evaluating the performance of ML models ranging from 1–10 predictors to validate the robustness of these biomarkers. It is important to keep in mind the relatively small sample size, which may affect the generalizability of our findings. This constraint underscores the need for further studies with larger cohorts to validate and expand upon our results.

Our findings were particularly significant in the evaluation of predictive performance across various ML models. When narrowed down to an optimized set of eight specific genes from our top 10, we observed enhanced prediction accuracies. The accuracy levels achieved by different models were as follows: LOR at 99.67 ± 2.36%, RF at 95.67 ± 12.05%, K-NN at 66.33 ± 29.45%, SVM at 88.00 ± 19.06%, and NN at an impressive 99 ± 4.00%. These findings emphasize both the robustness and reliability of our chosen genetic markers and underscore their clinical effectiveness in distinguishing between ‘SUICIDE’ and ‘NON-SUICIDE’ BD patients. The trained models were used finally to assess 7 BD patients, previously un-diagnosed but identified as high-risk due to factors like family history of suicide and personal attempts. Using a majority decision from the classifiers, the analysis indicated that two patients, UOHL027 and UOHL031, are at a low risk of suicide, while the remaining five are at high risk. Interestingly, UOHL027 has never attempted suicide, nor did he have any family history of suicide. Patient UOHL031 attempted suicide once and had one family member who died by suicide. Most high-risk patients either had multiple personal attempts or multiple family members with a history of suicide (except for patient UOHL030, who had one only one family member). This implies that our findings could contribute to the development of targeted therapeutic strategies and improve preventive measures in BD management. These aspects are crucial for enhancing patient outcomes and could potentially lead to personalized treatment approaches based on genetic predispositions. Our research has shown that the use of iPSC models can significantly advance the study of BD, potentially leading to personalized treatment strategies and improved preventive measures. These models allow for an in-depth analysis of BD at the cellular and molecular levels using neuronal cell types derived from patients. Findings from these models have highlighted mitochondrial dysfunction, neuron hyper-excitability, and disrupted calcium signaling, as well as distinct electrophysiological differences in the hippocampal neurons of patients responsive and non-responsive to lithium treatment [144,145,146,147,–148]. Building upon these insights, we further conducted a comparative analysis of mitochondrial and transporter gene expression, which revealed significant dysregulation of both mitochondrial and transporter genes in our study (Supplementary Figs. 6 and 7). These results underscore the multifaceted molecular alterations underlying BD and suicidality and highlight the relevance of mitochondrial and transporter pathways as potential contributors to disease pathophysiology and therapeutic targets. This research has contributed to the development of computational models that predict lithium treatment outcomes [53, 72]. Given the high risk of suicide among patients with BD, there is a strong need to understand the neurobiological and genetic factors of suicide, an iPSC model may shed more light on the related mechanism [149]. We believe the integration of ML with transcriptomic biomarkers is a promising and scalable approach, with the potential to extend suicide risk prediction beyond clinical populations in future research.

In conclusion, our study provides new insights into the genetic underpinnings of suicide in BD, elucidating complex interactions among genes involved in primary immunodeficiency, cellular signalling pathways such as cell adhesion and ion channels, and their links to other conditions like cardiovascular disorders. Through rigorous analysis of RNA sequencing data, we identified distinct gene expression signatures associated with low and high risks of suicide in BD patients and constructed high-accuracy classifiers, highlighting potential biomarkers that could lead to more precise, and accessible diagnostic methods. While this study offers important advances, the relatively small cohort size may limit the statistical power and generalizability of the findings. Additionally, although we were able to include detailed data on many clinical features, comprehensive information on variables such as rapid cycling, certain comorbidities (e.g., substance use, trauma exposure), and other potential confounders was not available for all participants. Despite these constraints, our study leverages a well-characterized cohort with detailed molecular profiling and addresses a critical gap in the literature. Future research using larger and more diverse samples with expanded clinical data will be essential to build on these findings and further elucidate the biological mechanisms underlying suicide risk in BD.

Correspondence and requests for materials should be addressed to Prof. Shani Stern. Additional information about the samples will be available in the GitHub file at this link (https://github.com/omveersharmanet/Predicting-Suicide-Risk-in-Bipolar-Disorder patients.git).

References

Author notes

1. These authors contributed equally: Omveer Sharma, Ritu Nayak.

Authors and Affiliations

1. Sagol Department of Neurobiology, University of Haifa, Haifa, Israel

   Omveer Sharma, Ritu Nayak, Liron Mizrahi, Wote Amelo Rike, Ashwani Choudhary, Hagit Sadis, Yara Hussein, Idan Rosh, Utkarsh Tripathi, Aviram Shemen, Yam Stern & Shani Stern

2. Department of Biomedical Sciences, University of Cagliari, Cagliari, Italy

   Alessio Squassina

3. Department of Psychiatry, Dalhousie University, Halifax, NS, Canada

   Martin Alda

4. National Institute of Mental Health, Klecany, Czech Republic

   Martin Alda

Contributions

Omveer Sharma (OS) and Ritu Nayak (RN) (co-first authors) contributed equally to this work. OS performed the data analysis, conducted machine learning predictor analysis, and drafted and compiled the manuscript. RN grew the cells, performed RNA extraction experiments, and drafted and compiled the manuscript. Liron Mizrahi (LM) assisted with Machine learning code. Wote Amelo Rike (WAR) and Ashwani Choudhary (AC) assisted with cell culture. Hagit Sadis (HS) helped in preparation of figures. Yara Hussein (YH) and Aviram Shemen (AS) assisted with data analysis. Idan Rosh (IR) provided cell care. Utkarsh Tripathi (UT) and Yam Stern (YS) reviewed the manuscript. Alessio Squassina (AS) provided the clinical samples for the cohort study. Martin Alda (MA) contributed to the study design and provided resources. Shani Stern (SS) (corresponding author) conceptualized and supervised the study, analyzed the data, and drafted the manuscript. All authors reviewed the manuscript.

Corresponding authors

Correspondence to Alessio Squassina, Martin Alda or Shani Stern.

Competing interests

The authors declare no competing interests

Ethics approval and consent to participate

All methods were performed in accordance with the relevant guidelines and regulations. All participants provided informed consent before participating in the study. Approval for the study was obtained from the Research Ethics Board of the University of Cagliari, Italy.

Open Access This article is licensed under a Creative Commons Attribution-NonCommercial-NoDerivatives 4.0 International License, which permits any non-commercial use, sharing, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if you modified the licensed material. You do not have permission under this licence to share adapted material derived from this article or parts of it. The images or other third party material in this article are included in the article’s Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article’s Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by-nc-nd/4.0/.

Reprints and permissions