Harnessing agent-based frameworks in CellAgentChat to unravel cell–cell interactions from single-cell and spatial transcriptomics
- 1School of Computer Science, McGill University, Montreal, Quebec H3A 2A7, Canada;
- 2Meakins-Christie Laboratories, Translational Research in Respiratory Diseases Program, Research Institute of the McGill University Health Centre, Montreal, Quebec H4A 3J1, Canada;
- 3Mila - Quebec Artificial Intelligence Institute, Montreal, Quebec H2S 3H1, Canada
Abstract
Understanding cell–cell interactions (CCIs) is essential yet challenging owing to the inherent intricacy and diversity of cellular dynamics. Existing approaches often analyze global patterns of CCIs using statistical frameworks, missing the nuances of individual cell behavior owing to their focus on aggregate data. This makes them insensitive in complex environments where the detailed dynamics of cell interactions matter. We introduce CellAgentChat, an agent-based model (ABM) designed to decipher CCIs from single-cell RNA sequencing and spatial transcriptomics data. This approach models biological systems as collections of autonomous agents governed by biologically inspired principles and rules. Validated across eight diverse single-cell data sets, CellAgentChat demonstrates its effectiveness in detecting intricate signaling events across different cell populations. Moreover, CellAgentChat offers the ability to generate animated visualizations of single-cell interactions and provides flexibility in modifying agent behavior rules, facilitating thorough exploration of both close and distant cellular communications. Furthermore, CellAgentChat leverages ABM features to enable intuitive in silico perturbations via agent rule modifications, facilitating the development of novel intervention strategies. This ABM method unlocks an in-depth understanding of cellular signaling interactions across various biological contexts, thereby enhancing in silico studies for cellular communication–based therapies.
Cells, as the fundamental units of most life forms, are dynamic and continuously adapt their states in response to changes in their microenvironment (Hu and Polyak 2008; Shao et al. 2020). In multicellular organisms, cellular activities rely on the coordination of cell–cell interactions (CCIs) across different cell types and tissues (Rouault and Hakim 2012; Bonnans et al. 2014; Zhou et al. 2018; Armingol et al. 2021). CCIs are critical in numerous biological processes, including cell differentiation (Vento-Tormo et al. 2018), disease pathogenesis (Stubbington et al. 2017), and immune responses (Giladi et al. 2020). By studying these interactions, we can uncover mechanisms behind biological processes, informing research on disorders and aiding in the discovery of diagnostic and therapeutic strategies (Giladi et al. 2020).
Inferring CCIs in biological processes is challenging owing to cellular heterogeneity in most biological systems (Yuan et al. 2017). Despite advancements in single-cell sequencing technologies like single-cell RNA sequencing (scRNA-seq) and spatial transcriptomics measurements (Eng et al. 2019; Rodriques et al. 2019; Xia et al. 2019; Marx 2021; Chen et al. 2022; Janesick et al. 2023), existing methods for CCI inference face significant drawbacks. Tools, such as CellPhoneDB (Efremova et al. 2020), CellChat (Jin et al. 2021), Connectome (Raredon et al. 2022), SingleCellSignalR (Cabello-Aguilar et al. 2020), and NATMI (Hou et al. 2020), rely on mean expression values from single-cell clusters, focusing on cell population levels and losing single-cell resolution. Scriabin (Wilk et al. 2024) analyzes LR interactions at the single-cell level but primarily relies on unimodal scRNA-seq data, limiting CCI inference ability owing to lack of spatial information. Because CCIs in tissue environments are heavily influenced by spatial structures and can occur at both long and short ranges (Asp et al. 2020), integrating spatial information is essential for studying endocrine and paracrine signaling. To address this, recent methods such as COMMOT (Cang et al. 2023), NICHES (Raredon et al. 2023), and CellChat v2 (Jin et al. 2025) have been developed for inferring CCIs from single-cell spatial measurements. However, no method effectively models both long-range and short-range interactions. Finally, CCIs play crucial roles in the pathogenesis of many diseases (Stubbington et al. 2017), and the analysis of these interactions holds the potential for developing effective therapeutics. Yet, existing methods like COMMOT and scSeqComm (Baruzzo et al. 2022) do not support in silico perturbation to utilize these impacts for exploring potential treatments, such as receptor antagonists.
Agent-based models (ABMs) simulate the actions and interactions of autonomous agents to understand complex system behaviors (Bonabeau 2002). Agents in ABMs represent entities with specific characteristics, behaviors, and interaction rules. ABMs have several advantages for modeling biological systems. First, ABMs capture emergent phenomena arising from collective agent interactions, revealing insights into cellular community function (Bonabeau 2002; Yu and Bagheri 2020). Second, they model heterogeneous cells as agents with distinct statuses, capturing variations in cell types, states, and behaviors (Bonabeau 2002; Yu and Bagheri 2020). Lastly, ABMs facilitate the straightforward and effective adjustment of agent rules, enabling perturbations in agent behavior.
To bridge the aforementioned gaps, we introduce CellAgentChat, a comprehensive ABM framework for inferring and visualizing CCIs from scRNA-seq and spatial transcriptomics data. Leveraging the advantages of ABMs, CellAgentChat offers versatile functionalities for CCI inference and visualization. By adopting an agent-based perspective, CellAgentChat integrates spatial and biological data, providing a more precise and comprehensive understanding of cellular interaction dynamics at both individual cell and system levels. Moreover, CellAgentChat enables in silico perturbation analysis of the downstream effects of cellular interactions on gene expression, facilitating the identification of potential therapeutic receptor targets. Notably, the adaptive framework of CellAgentChat allows for diverse functionalities by adjusting the behavior rules associated with cell agents. This feature enables the identification of long-range and short-range LR interactions. Finally, CellAgentChat can analyze temporal dynamics, shedding light on how CCIs influence cellular state and behavior over short time frames. Our evaluation of CellAgentChat, using several publicly available scRNA-seq and spatial transcriptomics data sets, reveal superior CCI inference compared to other benchmarked methods. Together, we introduce a new framework based on ABMs for future advancements in the study of cellular interactions implicated in vital biological processes.
Results
Overview of CellAgentChat
CellAgentChat constitutes a comprehensive ABM framework integrating gene expression data and existing knowledge of signaling LR interactions to compute the probabilities of CCI. Utilizing the principles of ABM, we characterize each cell agent through various cellular attributes (states), including cell identities (e.g., cell type), gene expression profiles, LR universe, and spatial coordinates (Fig. 1A). CCIs are measured by the quantity of ligands emitted by sender cells and subsequently absorbed by receiver cells. This process hinges upon three interrelated components: ligand diffusion rate (γl), receptor receiving rate (αr), and receptor conversion rate (βr) (Fig. 1B–D).
Overview of CellAgentChat. (A) CellAgentChat requires the following agent states for each cell: gene expression from scRNA-seq, cell type/cluster information, and LR database. Spatial coordinates (from spatial transcriptomics data) are optional. (B) The receptor receiving rate is the probability that an interaction is received by a receptor, calculated by non-parametric min-max normalization of receptor expression. (C) Our deep learning model inputs the gene expression of ligands and receptors, predicting the expression of all genes. The model incorporates prior knowledge of TF interactions. Receptor blocking is simulated by permuting receptor features and scaling down expression. The receptor conversion βr of the receptor is estimated from the total change of the predicted expression after the perturbation. (D) LR interaction scores depend on ligand diffusion rate (γl), receptor receiving rate (αr), receptor conversion rate (βr), and optionally spatial distance (d) between cells. Hyperparameters τ (spatial freedom, default = 2) and δ (ligand decay rate, default = 1) control interaction range. Using the LR score, we can calculate the interaction score between two clusters/cell types and the cell receiving score (CRS) for a cell. After a permutation test, significant interactions are identified and used to calculate the cell type pair score (CTPS) between two clusters. (E) CellAgentChat enables agent-based animation to visualize cell communication, identify significant LR pairs, and simulate receptor blocking to find the most perturbed downstream genes for potential therapeutic drug discovery.
The initial step involves the estimation of the receptor receiving rate (αr), which signifies the fraction of secreted ligands received by the receptor on the receiver cell (Fig. 1B). We then leverage a neural network to evaluate the receptor conversion rate (βr), which indicates how strongly an interaction influences downstream gene expression (Fig. 1C). Considering the possibility that a similar interaction strength between different LR pairs can yield diverse impacts on downstream gene expression dynamics (Armingol et al. 2021), this step is essential. Next, we compute the ligand diffusion rate (γl), which quantifies the ligands secreted by the sender cell that reach the target receiver cell (Fig. 1D; Eq. 1 in Methods). This rate is a function of both ligand expression, xl, intercellular distance, d, and two parameters: τ and δ. τ represents the degrees of spatial freedom concerning the single-cell data, typically set to two for spatial transcriptomics data that is derived from two-dimensional slices. δ signifies the decay rate of ligand diffusion, which prioritizes interactions over long or short distances. CellAgentChat quantifies interactions at multiple levels: between individual cells using LR scores, per-cell receiving strength via the cell receiving score, and between cell type populations. Together, CellAgentChat offers a suite of downstream functionalities (Fig. 1E), several of which are absent in existing methods.
Benchmarking CellAgentChat's CCI inference across diverse application scenarios
We undertook a comprehensive benchmarking of CellAgentChat against nine prominent CCI inference methods with the help of the LIANA+ benchmark (Dimitrov et al. 2024): CellPhoneDB (v5) (Efremova et al. 2020), CellChat (v2) (Jin et al. 2021, 2025), NICHES (Raredon et al. 2023), COMMOT (Cang et al. 2023), Scriabin (Wilk et al. 2024), Connectome (Raredon et al. 2022), SingleCellSignalR (Cabello-Aguilar et al. 2020), NATMI (Hou et al. 2020), and scSeqComm (Baruzzo et al. 2022). We evaluated the performance of all these methods across five different spatial transcriptomics data sets at both single-cell (or subcellular) and spot-level resolution: Stereo-seq mouse developing hippocampus at three developmental time points (E12.5, E14.5, E16.5) (Chen et al. 2022), 10x Genomics Visium human squamous cell carcinoma (Ji et al. 2020), 10x Genomics Visium pancreatic ductal adenocarcinoma (Moncada et al. 2020), SeqFISH+ mouse somatosensory cortex (Eng et al. 2019), and Xenium human breast cancer (Janesick et al. 2023). To deliver a comprehensive benchmarking, we compared all methods on these data sets from several key perspectives.
First, we compared the agreement of the cell–cell communication networks between each benchmarked method. We began by benchmarking across methods using only scRNA-seq data to validate the robustness of our CCI calculation framework before integrating spatial data. We assessed the cell–cell communication networks generated by CellAgentChat and compared them with the outcomes of other methods. Because of the absence of a definitive ground-truth, we adopted a consensus approach. This involved normalizing the interaction matrices computed by each method and averaging them to create an ensemble interaction matrix, serving as our new consensus reference. Subsequently, we computed the Pearson's correlation between the ensemble interaction matrix and the inferred communication patterns of each method for comparison. Across all data sets, the results from CellAgentChat demonstrated superior performance, achieving a higher median Pearson's correlation compared to all methods (Fig. 2A, top). Notably, CellAgentChat showed significantly better results than CellPhoneDB, NICHES, SingleCellSignalR, NATMI, and scSeqComm, further underscoring its consistency and accuracy in inferring communication networks (Fig. 2A, top). On examining individual data sets, CellAgentChat consistently outperformed CellPhoneDB, scSeqComm, and Scriabin in all cases (Wilcoxon one-sided test; P < 0.05) (Supplemental Figs. 1–3). Additionally, CellAgentChat and CellChat consistently showed high performance, but CellAgentChat exhibited less variability across all evaluated data sets (Fig. 2A; Supplemental Figs. 1–3). Upon incorporating spatially resolved transcriptomics data, CellAgentChat continued to perform with a high degree of consistency among the state-of-the-art methods. Specifically, CellAgentChat maintained superior performance compared to CellPhoneDB and also outperformed CellChat (Wilcoxon one-sided test; P < 0.05) (Fig. 2A, bottom; Supplemental Figs. 1, 4, 5). NICHES had a slightly higher median score (Wilcoxon one-sided test; P = 0.76), and COMMOT showed a slightly lower median (Wilcoxon one-sided test; P = 0.53), but neither difference reached statistical significance. Thus, CellAgentChat performed comparably to these methods, consistently exhibiting lower variability (9.9% smaller interquartile range compared to NICHES and 10.9% compared to COMMOT) and reinforcing its position among the top-performing methods.
Benchmarking of CellAgentChat with existing state-of-the-art methods. (A) Box plot comparing the Pearson's correlation between inferred cell communication by individual methods and the ensemble of all methods, without (top) and with (bottom) spatial information across all data sets. P-values calculated using Wilcoxon one-sided test: (*) P < 0.05; (**) P < 0.01; (***) P < 0.001. (B) Box plot comparing the Pearson's correlation between the inferred cell communication for each method and the inverse cell-distance matrix (top), across data sets. P-values calculated using Wilcoxon one-sided test: (*) P < 0.05; (**) P < 0.01; (***) P < 0.001. Bottom plot in panel shows the average ranking (lower is better) of each method's correlation across data sets. (C) The Reactome pathway enrichment analysis, performed on ligands and receptors sourced from the top 100 LR pairs inferred by each method across various data sets using spatial transcriptomics data. P-values calculated using Wilcoxon one-sided test: (*) P < 0.05; (**) P < 0.01; (***) P < 0.001.
Secondly, we examined the consistency of CellAgentChat's predictions with spatial information, focusing on the spatial proximity of stronger interactions. Specifically, we measured the Pearson's correlation between the probability of interaction and the inverse distance between cells, grounding this analysis on the established understanding that cells are more likely to interact when they are in closer proximity (Armingol et al. 2021). We excluded Connectome, SingleCellSignalR, NATMI, scSeqComm, and Scriabin for this comparison as they lack spatial data support. Indeed, cells predicted by CellAgentChat as highly interacting were closer in proximity compared to predictions from other methods (Fig. 2B). We demonstrated that our spatially informed predictions outperformed those of CellPhoneDB and CellChat (Wilcoxon one-sided test; P < 0.05) (Fig. 2B), showing superior performance across every data set compared to CellPhoneDB and in all but one data set compared to CellChat (Supplemental Fig. 6). Furthermore, we show comparable performance to NICHES and COMMOT across data sets (Fig. 2B; Supplemental Fig. 6). Overall, CellAgentChat and COMMOT had the best average rank (lower is better) across all spatial data sets (Fig. 2B, bottom).
Benchmarking CellAgentChat's CCI predictions via pathway enrichment analysis and ligand–receptor prediction agreement
After validating the accuracy of CellAgentChat in predicting CCIs through consensus and spatial proximity metrics, we conducted Reactome pathway enrichment analysis to further confirm the biological relevance of the inferred LR pairs. This step serves as an additional layer of validation. We conducted computational validations comparing the enrichment P-values for relevant pathways using ligands and receptors from the top 100 significant LR pairs from each method.
The pathway enrichment analysis conducted for the results using only scRNA-seq data underscored CellAgentChat's strength over other methods. CellAgentChat consistently outperformed Scriabin across all data sets (Supplemental Fig. 7). Compared to CellPhoneDB and SingleCellSignalR, CellAgentChat achieved superior performance in six of the seven data sets and comparable results in the remaining one. Additionally, it outperformed scSeqComm and NATMI in the majority of data sets, while showing comparable performance in the rest (Supplemental Fig. 7). Upon integrating spatial data, pathway enrichment analysis continued to highlight CellAgentChat's superiority over CellPhoneDB and CellChat (Fig. 2C). CellAgentChat continues to perform better than CellPhoneDB in all data sets (Fig. 2C; Supplemental Fig. 7B). CellAgentChat also demonstrated superior performance to CellChat, NICHES, and COMMOT in the majority of data sets (Fig. 2C; Supplemental Fig. 7B). Overall, pathway enrichment analysis supported the CCI inference precision of CellAgentChat over other benchmarked methods.
Finally, we calculated the overlap in LR pairs among the top 100 identified LR pairs across the nine methods. We observed that the overlap between CellAgentChat and the other methods, particularly CellChat, NICHES, and COMMOT, was high (Supplemental Figs. 8, 9). This further confirms CellAgentChat's effectiveness, as its LR signaling predictions align with state-of-the-art methods.
Building on the findings from the preceding benchmarking sections, which underscore CellAgentChat's superior, or at least comparable, performance compared to other methods, it is crucial to recognize that beyond these achievements, CellAgentChat also introduces distinct functionalities not available in other tools (Table 1). Unique to CellAgentChat is its integrated ABM animation platform, which facilitates the visualization of individual-level cellular interactions and supports dynamic simulations. Additionally, through easy and flexible adjustment of agent behavior, CellAgentChat empowers the use of in silico receptor blocking to explore potential interventions effectively and serves as a valuable tool for identifying both long- and short-range interactions. These functionalities will be discussed in more detail in the following sections, showcasing how they contribute to the broader utility and potential of CellAgentChat.
Summary of properties and functionalities for all benchmarked CCI inference methods
CellAgentChat identifies signaling pathways using spatial transcriptomics data at single-cell and multi-cell/spot resolutions
To demonstrate CellAgentChat's accuracy in predicting CCIs, we applied it to single-cell spatial transcriptomics data from the Stereo-seq platform during mouse hippocampus development (Chen et al. 2022). We focused on identifying key signaling pathways involved in the developmental lineage between cells at embryonic day E12.5, E14.5, and E16.5 (Fig. 3A).
Analysis of cell communications in developing mouse hippocampus and squamous cell carcinoma (SCC) using CellAgentChat. (A) UMAP of scRNA-seq data across time points for the developing mouse data set. (B) Heat map displays inferred cell communication network at E12.5 in the developing mouse. (C) Dot plot of the top 25 LR pairs identified at E12.5 in the developing mouse. (D) Heat map displays inferred cell communication network in the SCC. (E) Dot plot of the top 25 LR pairs identified in SCC.
CellAgentChat identified 1531 significant interactions involving 579 LR pairs at E12.5 (Fig. 3B,C). Microglial (Micro), fibroblasts (Fibro), erythrocytes (Ery), endothelial cells (Endo), and glioblasts (GlioB) had the highest number of interactions (Fig. 3B), aligning with their known role in neurogenesis (Tong and Vidyadaran 2016; Li et al. 2022; Vogenstahl et al. 2022). Key LR pairs identified included the Apoe signaling pathway (Fig. 3C), consistent with previous findings on the regulatory role of Apoe in neurogenesis and regulating neurodegenerative diseases (Tensaouti et al. 2018; Yin et al. 2023). Additionally, we identified significant collagen signaling (Col1a1, Col1a2, and Col8a1), promoting hippocampal neurogenesis in mice (Kakoi et al. 2012; Swonger et al. 2016). To systematically assess the reliability of the LR interactions predicted by CellAgentChat, we employed Reactome pathway enrichment analysis. Our analysis revealed significant enrichment of pathways related to neurogenesis, including axon guidance, nervous system development, developmental biology, and extracellular matrix (ECM) organization (Cope and Gould 2019; Fromme and Zigrino 2022), with CellAgentChat outperforming existing methods (Fig. 2C). Further analyses at E14.5 and E16.5 yielded consistent results (Supplemental Results; Supplemental Figs. 10, 11).
CellAgentChat seamlessly integrates multiple spatial sections into a single analysis, as demonstrated with three combined slices at E16.5 (Supplemental Fig. 10B). This reveals consistent communication patterns across all slices (Supplemental Fig. 12A–C), demonstrating the tool's robustness and efficacy in capturing coherent signaling networks (Supplemental Fig. 12D). To further showcase broad applicability, we also applied CellAgentChat without spatial information, obtaining consistent results (Supplemental Fig. 11).
To showcase the broad applicability of CellAgentChat, we applied it to spot-level Visium spatial transcriptomics data from squamous cell carcinomas (SCC) (Ji et al. 2020), with a total of 13 cell populations, including normal keratinocytes (NKer), tumor keratinocyte (TKer), and myeloid dendritic cell (MDC) (Fig 3D). Treating each spot as a cell, CellAgentChat uncovered key LR pairs driving SCC progression, notably PKM-CD44, MMP1-CD44, and HLA-B-CANX, consistent with prior studies (Fig. 3E; Pontén et al. 2008; Atasoy et al. 2009; Ludwig et al. 2019; Schaafsma et al. 2021; Boschert et al. 2022; Zhang et al. 2022). Reactome pathway analysis revealed significant enrichment of pathways relevant to SCC, including ECM organization (Fromme and Zigrino 2022; Parker et al. 2022), immune system (Ahrends and Borst 2018), MET signaling (Szturz et al. 2017; Centuori and Bauman 2022), and laminin interactions (Fig. 2C; Ramovs et al. 2017; Meireles Da Costa et al. 2021). CellAgentChat outperformed other methods in identifying these pathways, highlighting its superior ability to predict biologically relevant interactions in cancer (Fig. 2C).
To further demonstrate the broad applicability of our model, we also applied CellAgentChat to diverse single-cell spatial data sets across different platforms, including Xenium breast cancer, SeqFISH+ mouse somatosensory cortex, and human pancreatic ductal adenocarcinoma (PDAC) (Supplemental Figs. 13, 14), thereby confirming its compatibility and functionality across diverse spatial technologies.
CellAgentChat enables in silico receptor block perturbation to identify potential therapeutics
CellAgentChat facilitates in silico receptor blocking to simulate the interruption of specific signaling pathways. By computationally blocking receptors, researchers can observe and predict the impact on downstream genes, which may be associated with disease processes. This approach allows for the identification of potential therapeutic receptor targets without the need for physical experiments. Such in silico perturbations offer a methodical approach to therapeutic discovery, providing valuable insights into the molecular mechanisms of diseases and accelerating the identification of potential treatments.
We harnessed CellAgentChat to build upon previous research on breast cancer exploring the Xenium human breast cancer data set (Janesick et al. 2023), aiming to identify significant receptors related to the disease and cultivate innovative therapeutic approaches. The first step of our process entailed using CellAgentChat to pinpoint 13 key receptors involved in the top 25 LR interactions for the cellular communications in the data set (Supplemental Figs. 13D, 15A). We then simulated the blocking of these receptors individually, akin to the actions of receptor inhibitor drugs. We deployed a feedforward neural network that was designed to calculate each receptor's conversion rate and was trained to predict the expression of all genes using the ligand and receptor expression as inputs. The network connections were constrained to the known LR interactions, LR transcription factor (TF), and TF genes along the forward pass, respectively (Fig. 4A). Using the trained network, we predicted the impact of receptor blocking on downstream gene expression and identified the most perturbed downstream genes for each in silico blocked receptor. We subsequently assessed these perturbed genes using DisGeNET (Piñero et al. 2015) enrichment analysis to evaluate the potential therapeutic implications of blocking specific receptors (Fig. 4A). Our analysis with CellAgentChat indicated that blocking nine of the 13 identified top interacting receptors in the breast cancer data set—specifically EGFR, PDCD1, PECAM1, SDC4, CXCR4, AVPR1A, CTLA4, CD80, and CD93 (hereafter referred to as candidate receptors)—significantly disrupts genes associated with breast cancer (Fig. 4B).
CellAgentChat enables in silico receptor blocking to identify potential receptor antagonists for therapeutics. (A) Schema of the in silico perturbation analysis using our neural network to block receptors and identify the most perturbed downstream genes. Blue bars show unperturbed predicted expression; red bars show perturbed predictions. (B) Breast cancer gene set enrichment (−log10 adjusted P-values; binomial test, FDR corrected) for significantly perturbed genes after blocking 13 selected receptors. Receptors highlighted in red are significantly enriched. (C–E) Matrix plots showing mean expression changes per cell type for EGFR, PDCD1 (PD-1), and CTLA4, respectively (left). Top 10 downstream genes with highest expression change after blocking EGFR, PDCD1 (PD-1), and CTLA4, respectively; red highlights indicate breast cancer-associated genes (middle). Survival analysis comparing high versus low receptor expression groups for EGFR, PDCD1 (PD-1), and CTLA4, respectively (right).
In an ablation study to highlight the benefits of integrating biological priors, we tested a fully connected feedforward network without the LR TF and TF target prior constraints. The removal of these biological constraints resulted in a noticeable decrease in specificity: four receptors previously identified as nonsignificant were now marked as significant (Supplemental Fig. 15B).
CellAgentChat's analysis reveals that specific cell types, including tumor cells, stromal cells, endothelial cells, and macrophages, are highly engaged in interactions with the candidate receptors (Supplemental Fig. 16). By simulating the blocking of these receptors, CellAgentChat predicts a significant impact on these cellular interactions, demonstrating the tool's potential to identify effective therapeutic strategies.
Although studies have associated all of the candidate receptors examined in this study with breast cancer, the candidate receptors we identified—EGFR, PDCD1 (also known as PD-1), and CTLA4—appear to play a particularly crucial role in breast cancer (Fig. 4C–E; Masuda et al. 2012; Ralser et al. 2021; Ghahremani Dehbokri et al. 2023). Furthermore, therapeutic treatments targeting the three receptors in breast cancer have been previously documented. Both lapatinib and pembrolizumab are FDA-approved drugs designed to specifically target and disrupt pathways involving EGFR and PDCD1, respectively. Ipilimumab and tremelimumab are approved by the FDA to suppress the role of CTLA4 in melanoma (Ghahremani Dehbokri et al. 2023). Although ipilimumab and tremelimumab are not FDA-approved for breast cancer treatment, studies have shown their positive impact in this context (Ghahremani Dehbokri et al. 2023). This corroborates the potential application of CellAgentChat to identify receptor-mediated therapeutics. Moreover, blocking the EGFR, PDCD1, and CTLA4 receptors in silico perturbs the genes associated with breast cancer (from DisGeNET) (Fig. 4C–E). Our analysis revealed that blocking these receptors perturbs several disease genes, including TPD52, CDH1, and EPCAM, among others, which have been previously linked to breast cancer (Soysal et al. 2013; Corso et al. 2018; Ren et al. 2021). Numerous perturbed genes notably influenced immune cell types like mast cells and macrophages (Fig. 4C–E, left), which play a crucial role in the tumor microenvironment and immune response against cancer cells. All other breast cancer disease target genes perturbed are highlighted in red in the figure (Fig. 4C–E; Supplemental Figs. 15C–E, 17A–C).
To validate the accuracy of the candidate drug targets captured by CellAgentChat, we conducted a survival analysis of breast cancer patients over time, comparing populations with high and low expression levels of each candidate therapeutic receptor. The analysis revealed statistically significant differences in overall survival probabilities over 150 months for patients with varying levels of expression of EGFR, PDCD1, and CTLA4 (Kaplan-Meier log-rank test; P = 1.3 × 10−10, P = 3.4 × 10−8, P = 1.6 × 10−9, respectively) (Fig. 4C–E). Additionally, survival analysis of all other candidate receptors exhibited statistically significant differences in survival probabilities between the two patient groups (Supplemental Figs. 15C–E, 17A–C). These findings indicate that the identified candidate receptors would serve as suitable therapeutic targets for treating breast cancer, highlighting the powerful functionality of CellAgentChat's in silico receptor blocking to predict and validate potential therapeutic interventions.
To demonstrate the applicability of CellAgentChat's in silico perturbation on other spatial transcriptomics technology, we also demonstrate the in silico perturbation capabilities of our model on a spot-level resolution PDAC data set (Moncada et al. 2020), where we identified NOTCH3, PLAUR, ITGA11, and CD74 as potential therapeutic targets (see Supplemental Results; Supplemental Figs. 18, 19).
To validate the effectiveness of CellAgentChat in predicting the effects of receptor perturbation on disease target genes, we compared our approach with scGen (Lotfollahi et al. 2019). Using the breast cancer data set, we found that scGen did not identify any receptors whose perturbation significantly disrupted genes associated with breast cancer (Supplemental Fig. 20A). Additionally, receptor perturbation using CellAgentChat showed a higher specificity in affecting breast cancer-related genes compared to scGen, especially for well-known breast cancer receptors EGFR and ERBB2 (Supplemental Fig. 20B). This validation was also conducted for the PDAC data set (Supplemental Results; Supplemental Fig. 20C–E).
Analyzing perturbation effects in both diseased and normal tissues provides critical insights into underlying communication networks. CellAgentChat's ABM capabilities revealed three receptors—LRP1, CAV1, and CD93—that significantly disrupt PDAC-associated genes when blocked in tumor regions, while having no effect in normal tissue (see Supplemental Results; Supplemental Figs. 21–25).
CellAgentChat facilitates agent-based visualization of interactions associated with individual cells
The ABM approach used in CellAgentChat allows us to examine individual cell interactions rather than solely focusing on cell populations. We applied CellAgentChat to a spot-level human PDAC data set (Moncada et al. 2020), with six distinct cell populations: ductal (Duct), acinar (Acin1, Acin2, and Acin3), stromal (Strom) cells, and two classes of cancer (Can1 and Can2) (Supplemental Fig. 13B,E,F). To demonstrate CellAgentChat's capability in revealing insights into cell interactions, we identified the 20 most receiving spots (group A, red circles) based on the cell receiving score (Methods) and the least receiving spots (group B, black circles) within Can2 cells (green cells) (Fig. 5A).
CellAgentChat facilitates inference of cellular interactions associated with individual cells in the PDAC data set. (A) CellAgentChat's animation platform enables the detection of communication strength variations in individual cancer cells. The platform highlights the cell receiving score (Methods) of Can2 cells (colored green). The top 20 cells with the highest CRS (group A) are marked by red circles and the bottom 20 cells with the lowest CRS (group B) are marked by black circles. (B) Expression analysis of receptors in the top 25 LR interactions shows significantly higher expression in group A compared to group B (Wilcoxon test one-sided; P < 0.05). (C) Mean expression of differentially expressed (DE) genes upregulated for group A spots relative to rest of the spots (top) and upregulated for group B cells relative to rest of the spots (bottom). (D) Heat maps display the DE genes for group A (left) and group B (right). The DE genes of group A cells are mainly expressed by group A and not by group B (Wilcoxon test one-sided; P < 0.001). Both groups express the DE genes of group B but are more significantly expressed by group A spots (Wilcoxon test one-sided; P < 0.05). (E) The results of the Reactome pathway analysis using group A DE genes (left) and group B DE genes (right). Red highlights indicate shared pathways.
The spots with higher CRS (more interactive) generally reside in the outer layer of the tissue (invasive front), consistent with existing studies (Friedl and Alexander 2011). These two cell groups of distinct interacting patterns within the same cluster also demonstrated significant transcriptomic differences. Group A showed higher mean expression levels of receptors involved in the top 25 LR pairs than group B (Fig. 5B). Additionally, receptors such as LRP1, SDC1, ITGA3, NOTCH3, and CD74 displayed notable differences in mean expression between the two groups, with group A again exhibiting higher expression (Fig. 5B).
Next, we systematically explored the differentially expressed (DE) genes of these two groups by comparing the spots within the group with the rest of the spots in the data set. The differential genes associated with group A show significantly higher expression in group A cells compared with group B (Wilcoxon one-sided test; P = 1.22 × 10−48) (Fig. 5C,D). In contrast, the expression of group B DE genes is relatively similar across groups (Fig. 5C,D). To identify the biological functions associated with these groups, we conducted Reactome pathway enrichment on the top 50 DE genes for both groups (Fig. 5E). Further, group A was associated with several pathways related to the extracellular matrix and collagen, known to be crucial in PDAC tumor progression (Perez et al. 2021). Altogether, CellAgentChat identifies complex interactions at the individual cell level, revealing diverse roles in disease progression, highlighting the tool's potential to uncover new insights in cellular communication.
CellAgentChat can be tuned to identify long- and short-range interactions
Our understanding of the specific distances at which LR pairs operate is evolving but remains incomplete. These interactions can occur closely between cells or span longer distances, influenced by factors like specific ligands/receptors, tissue environment, and physiological conditions. CellAgentChat's ABM captures detailed cell behaviors, including the diffusion of signaling ligands, enabling modeling of both short- and long-range interactions at the single-cell level by adjusting the delta (δ) factor to control the ligand diffusion rates (Methods).
Using CellAgentChat, we investigated the spatial distance tendencies of cellular communication within the mouse somatosensory cortex leveraging data from SeqFISH+ technology (Supplemental Fig. 14; Eng et al. 2019). To this end, we utilized the CellChat database, which categorizes LR pairs into “secreted signaling” for long-range interactions and “cell–cell contact” for short-range interactions based on known distance ranges.
To identify short-range interactions, we used a δ = 10 to constrain interactions to shorter distances and used δ = 0.1 to identify long-range interactions. To validate our results, we compared the top interacting LR pairs with CellChat's annotations. Under our short-range mode, the majority of our top 10, 15, 20. and 25 LR pairs were annotated as “cell–cell contact” (short-range) (Fig. 6A). Specifically, 80% of the top 15 LR pairs identified in the short-range mode fell into the “cell–cell contact” category (Fig. 6A). In the long-range mode, most of the top LR pairs identified were classified as “secreted signaling” (long-range), with 75% of the top 15 LR pairs falling into this category (Fig. 6A). Recently, Liu et al. (2022) have developed a method to identify long-range and short-range interactions using optimal transport (OT). Compared to their approach, CellAgentChat demonstrated superior accuracy in capturing short-range interactions and similar accuracy in identifying long-range interactions (Fig. 6B).
CellAgentChat identifies long-range and short-range interactions in the adult mouse cortex data set. (A) Bar plot illustrates the percentage of interactions found by CellAgentChat annotated by CellChat's database as either “cell–cell contact” (short-range) or “secreted signaling” (long-range) with a delta (δ) value set to 10 (top) and 0.1 (bottom). (B) Bar plot illustrates the percentage of interactions found by CellAgentChat annotated by CellChat's database as either “cell–cell contact” or “secreted signaling” in comparison to an optimal transport (OT)-based approach (Liu et al. 2022). (C) CellAgentChat's animation platform depicts the cell receiving score for a “cell–cell contact” Sema4a-Plxna4 ligand–receptor pair using the short-range mode (δ = 10). The color of the cells (circles) represents the cell cluster. This animation depicts the release of Sema4a (ligand) from cells in cluster 9 (represented by turquoise cells). In the animation, cell size corresponds to the cell receiving score. Cells closer to cluster 9 cells exhibit higher CRS (within red square) than cells that are further away (within blue square), which are smaller in size. (D) Reactome pathway analysis of the receptors involved in the top 25 LR interactions in the long-range mode (top) and short-range mode (bottom).
In our short-range mode, we detected LR pairs like Sema4a-Plxna4, known for its role in axon guidance in the mouse cortex, over short distances (Neufeld et al. 2016). We further employed our animation platform to visualize the CRS of individual cells for the Sema4a-Plxna4 LR pair in the short-range mode originating from cluster 9 cells (turquoise cells) (Fig. 6C). The results highlighted higher CRS among cells close in proximity to the cluster 9 cells represented by larger sizes (Fig. 6C, within the red box), compared to cells situated further away (Fig. 6C, within the blue box). This observation underscores the capabilities of CellAgentChat in identifying and modeling short-range interactions.
We conducted pathway enrichment comparisons of the LR interaction pairs predicted by CellAgentChat for both long- and short-range interactions (Fig. 6D). The genes associated with short-range interactions showed enrichment in pathways related to cell adhesion and cell–cell junctions. In contrast, genes associated with long-range interactions were enriched in numerous signaling pathways known for their wide regulatory range (Fig. 6D). This analysis underscores CellAgentChat's ability to distinguish between different types of cellular interactions and its broader utility in understanding how these interactions influence genomic processes.
To further validate the accuracy of CellAgentChat's predictions, we compared the molecular weights of ligands identified in the long-range mode (δ = 0.1) versus the short-range mode (δ = 10). Our analysis revealed statistically significant differences in the molecular weights of ligands across all top pair ranges (Supplemental Fig. 26). Specifically, ligands predicted in the long-range mode exhibited significantly lower molecular weights compared to those in the short-range mode. This finding aligns with expectations, as heavier ligands are less likely to diffuse over long distances and are thus more prevalent in short-range interaction (Li et al. 2009). In contrast, lighter ligands, which exhibit higher diffusivity, are more likely to mediate long-range interactions and are enriched in the long-range mode (Li et al. 2009). These results support the biological plausibility of our model's predictions and underscore the potential for integrating ligand molecular weight information into real-world data sets.
Ultimately, we can leverage our understanding of the identified short-range and long-range LR pairs to assign the ligand-specific
decay rates (δ), enabling precise control over ligand diffusion () and accurate modeling of both interaction types simultaneously.
CellAgentChat empowers dynamic simulations using agent-based modeling
CellAgentChat provides dynamic simulation capabilities to analyze how CCIs impact the cell states over short periods. At each time point (time step), CCIs are computed, followed by the prediction of gene expression levels for each cell using our trained neural network (Methods). Given uncertainties and various influencing factors affecting gene expression, directly substituting existing values with predicted ones is impractical. Instead, we opt for assessing the directional change of each gene—whether they increase or decrease—which proves to be more robust and easier to predict (Supplemental Fig. 27A; Ding et al. 2018). For each cell at every time step, genes undergoing significant changes are identified, and their expressions are adjusted in the simulation for subsequent steps (Methods; Supplemental Fig. 27A). This iterative process captures the dynamic modulation of gene expression by CCIs, providing insights into how cellular interactions influence genomic responses.
To validate the accuracy of our dynamic simulations, we determined the correlation between the pseudotime trajectory of each gene and its trajectory as established by the dynamic simulation of CellAgentChat across five time steps. Across all genes, we observed a correlation between the pseudotime trajectory and the ABM dynamic trajectory (Supplemental Fig. 27B). We found strong agreement between the pseudotime trajectory and the dynamic ABM trajectory for genes like ITGA6 (receptor) and MMP1 (ligand), which showed overall increasing trends over time (r2 = 0.93 and r2 = 0.97, respectively) (Supplemental Fig. 27C,D). Similarly, both trajectories indicated a decreasing expression trend for CALR (ligand) over time (r2 = 0.90) (Supplemental Fig. 27E). Moreover, CellAgentChat accurately identified the initial decrease of expression of CD44 (receptor) followed by a rapid increase (r2 = 0.89) (Supplemental Fig. 27F). These findings highlight the effectiveness of CellAgentChat in capturing dynamic cellular changes resulting from CCIs within a time frame comparable to the pseudotime represented by the cells. More validation into our dynamic simulations can be found in the Supplemental Results and Supplemental Figure 27.
Discussion
We introduced CellAgentChat, an agent-based model capable of elucidating cell–cell interactions from single-cell and spatial transcriptomics data. By modeling individual cell agents with simple behavior rules, CellAgentChat outperforms existing methods in reconstructing CCI networks, integrating spatial context, and identifying meaningful ligand–receptor (LR) interactions (Fig. 2). CellAgentChat, uncovered key LR pairs in neural development and cancer, such as Apoe and collagen signaling in mouse neurogenesis (Fig. 3), LAMC2-ITGA6 in squamous cell carcinoma (Fig. 3), and therapeutic receptor targets EGFR, PDCD1, and CTLA4 in breast cancer (Fig. 4). CellAgentChat also distinguishes long- and short-range interactions, identifying Sema4a-Plxna4 as a short-range axon-guidance signal in the mouse cortex (Fig. 6; Neufeld et al. 2016). These results underscore CellAgentChat's versatility and robustness in advancing our understanding of complex biological systems and its potential impact on therapeutic strategies in biomedical research.
In contrast to the existing methods for predicting CCI, CellAgentChat offers various distinctive features. Firstly, it offers a dynamic visualization of cellular interactions related to individual cell agents. This animation tool allows researchers to explore the intricate heterogeneity within cell populations, revealing insights at both the cluster and individual cell levels that were previously inaccessible (Fig. 5) Furthermore, the animation feature enables researchers to visualize and explore the dynamic impact of CCIs over short time frames. As cell agents interact, their states—such as gene expression profiles—continuously update in response to CCIs. This dynamic visualization provides insights into how cellular behaviors and interactions evolve over time, offering a powerful tool to study emergent phenomena and identify key patterns or critical transitions driven by CCIs in a visually intuitive manner. Secondly, CellAgentChat enables efficient in silico receptor blocking, facilitated by modulating the agent rules (for instance, receptor receiving rate), allowing for a systematic assessment of a receptor blocking's impact on downstream gene expression dynamics. The perturbation pipeline is robust, specifically targeting disease-related genes rather than random genes, ensuring relevance to the biological context. By investigating these genes, CellAgentChat can propose systematic strategies for potential therapeutic interventions, making it a valuable tool for exploring targeted treatments. Thirdly, CellAgentChat can be tuned to detect both long- and short-range interactions, capturing cellular communication across spatial scales. This flexibility allows researchers to study diverse biological processes, such as endocrine signaling at a distance or juxtracrine signaling through direct contact, providing insights into tissue development, immune responses, and disease progression.
In future work, we will explore model-level optimizations tailored specifically for spatial transcriptomics data. These improvements may enhance the model's ability to capture spatial relationships, leading to better biological interpretability in our downstream functionalities. We also plan to extend CellAgentChat by incorporating cell movement and cell death, enhancing the model's realism. Although technically challenging, another promising extension is representing ligands, receptors, and cofactor molecules as individual agents, enabling more comprehensive simulations of complex cellular interactions. This approach would offer deeper insights into cell-molecule dynamics and improve the efficiency of in silico perturbations like receptor blocking.
Because of low gene coverage in single-cell resolution spatial datasets (e.g., SeqFISH+, Xenium), spot-level technologies remain popular. CellAgentChat allows agents to represent individual cells or spots for flexibility. Although spot-level representation may overlook some CCIs within spots, it simplifies complex analyses. For higher resolution, tools like SpatialScope (Wan et al. 2023) or scSemiProfiler (Wang et al. 2024) can deconvolve spots into single cells, enabling CellAgentChat to infer interactions, perform in silico perturbations, analyze long/short-range interactions, and explore key LR pairs (see GitHub repository for details; Software availability).
In this study, we avoided imputation or denoising techniques, focusing instead on developing CellAgentChat to infer CCIs directly from raw normalized spatial transcriptomics data. However, imputation methods such as MAGIC (Van Dijk et al. 2018), can enhance gene expression quality. Interested researchers can preprocess data accordingly before applying CellAgentChat, highlighting our framework's flexibility and ease of use (see GitHub repository for guidance; Software availability).
CellAgentChat enhances the study of cellular interactions by providing a versatile ABM framework. It facilitates dynamic visualization, prioritization of interaction ranges and in silico receptor blocking, aiding in the identification of therapeutic targets. By integrating spatial information and accommodating various data resolutions, CellAgentChat offers a valuable tool for biomedical research by enhancing the understanding of disease mechanisms and developing therapeutic strategies through the lens of CCIs.
Methods
Ligand–receptor database
In this study, we used the ligand–receptor interactions from CellTalkDB (Shao et al. 2021), a curated database of almost 3400 human and over 2000 mouse literature-supported interactions, as the default LR interaction resource for our CellAgentChat method. Besides the default CellTalkDB LR interactions, CellAgentChat also allows for customized LR interactions that users obtain from other sources or manually curate. For our analysis of long-range and short-range interactions, we instead utilized the CellChat database, which classifies LR pairs into “secreted signaling” for long-range interactions and “cell–cell contact” for short-range interactions, based on established distance ranges.
Scoring cell–cell interactions with the agent-based model
To build the agent-based model, we need to define the states and behavior of each agent. Agents here represent the cells in the biosystem. The agent states describe the state and features associated with each cell, whereas the agent rules define how each agent behaves.
Cell agent state description
As we are primarily interested in studying the interactions between cells (or clusters), we describe each cell agent with characteristics and features related to cellular interactions. The remainder of this section provides the detailed cell agent state description. (1) Cell Agent Expression: The expression value of all genes, including ligands and receptors, is stored for each cell. The gene expression profile presents the quantitative description of the cell agent state, which is essential to infer the interaction behaviors. (2) Cell Agent Spatial Position: We use spatial transcriptomics data to record the position of each cell agent. All the cells are placed on the grid at their spatial coordinate (if provided). Otherwise, CellAgentChat randomly distributes all cells on the grid. (3) Cell Agent Batch: To utilize multiple samples (biological replicates or technical samples) concurrently, we designate the batch of each cell to confine interactions solely between cells of the same batch. (4) Cell Agent Ligand–Receptor Universe: Every cell agent has the same LR universe using interactions curated from CellTalkDB or user-customized LR interactions. Users can modify the LR universe by blocking a specific receptor or LR interaction (see “Agent-based dynamic simulations and animation of CCIs” in Methods). By doing so, they can remove all interactions associated with ligands or receptors (or specific LR interactions) from the universe of all cells. This characteristic of ABMs enables effective and flexible ligand/receptor perturbations.
Cell agent communication behavior rules
In this study, the cellular communications between agents (cells) depend on the states of the cells, which include the gene expression of the cell agent, the spatial position of the cell agent (optional), and the LR universe. To be more specific, the cellular interaction strengths (scores) between cell agents are quantified based on the set of defined rules as discussed below.
-
Ligand diffusion rate (
): The diffusion rate
estimates the rate of the diffusion of a ligand l in the microenvironment from sender cell i to receiver cell j. The diffusion rate depends on a few factors. First, the diffusion rate is proportional to the number of ligands a sender cell secretes into the environment. We can estimate this from the number of RNA molecules proportional to its gene expression.
Second, the diffusion rate is inversely proportional to the distance (of the receiver agent) from the sender cell agent. We calculate the Euclidean distance between two cells based on the spatial coordinates of all the cell agents. The hyperparameter, τ, is a dimensionality parameter representing the degrees of spatial freedom for spatial data. τ is typically set to two for spatial transcriptomics data derived from two-dimensional slices. However, with the advancement of technology and the increasing prevalence of three-dimensional slices (Rao et al. 2021), we can adjust the τ parameter accordingly (τ = 3). When spatial information is unavailable (e.g., in scRNA-seq data), τ is set to 0 so that the distance between all cells in the data will degenerate to 1 (dτ=0 = 1).
In the absence of precise information about the rates of ligand travel or the operational distance range of LR pairs, we employ the delta (δ) parameter as an approximation for the distance range of LR pairs and model the decay rate of ligands. A user-adjustable slider δ (between 0 and 10, default = 1) is available to control the decay rate of the diffusion along with the increase of distance, which allows for inferring long-range or short-range interactions. Keeping δ at the default value of 1 provides no bias toward long-range or short-range interactions. We modeled short-range interactions with a δ = 10 and long-range interactions with δ = 0.1. By investigating various δ values and the corresponding agent communication behaviors defined by these values, we can also determine a specific δ value for each LR pair. This approach avoids the assumption that all ligands share the same diffusion rate. The distance between sender cell i and receiver cell j is calculated as the Euclidean distance d(i, j). The expression of ligand l in sender cell i is denoted as
. Then, the diffusion rate
can be calculated as below:
(1)
To address the variability in distance scales across different spatial techniques and data sets, we normalize spatial distances within each data set to a uniform range of 0–100. This normalization ensures that communication strengths are comparable across data sets while preserving the relative spatial relationships. Furthermore, this normalization is crucial for visualizing cell interactions within our ABM framework (see “Agent-based dynamic simulations and animation of CCIs” in Methods), as maintaining consistent spatial scales enhances interpretability.
-
Receptor receiving rate (
): The receiving rate
associated with a receptor r represents the probability that an interaction is received by the receptor in cell j, and it is proportional to the expression of the receptor. The receptor receiving rate will be a value in the range (0,1) where 0 represents that the receptor cannot “receive” any ligands that reach it, and 1 denotes that reached ligands will all be “received” by the receptor. Here, we employed the non-parametric min-max normalization to calculate the receptor receiving rate for each receptor. Min-max normalization has been used in other methods for cell interaction inference (Tran et al. 2022; Peng et al. 2023, 2024; Feng et al. 2024), lending support to our choice. For each cell, the minimum expression value in the cell is subtracted from the expression value of each receptor in that cell and then divided by the difference between the maximum and minimum expression values in that cell (min-max normalization), enabling the identification of highly expressed receptors within each cell, thereby establishing a cell-specific metric.
-
Receptor transcriptome conversion rate (βr): CCIs play a crucial role in the regulation of gene expression, impacting the transcriptome and ultimately influencing cellular function and fate. In order to account for the effect of cellular interaction on downstream gene expression, we employed a neural network regressor to estimate the downstream gene expression impact associated with each receptor (here, we term it as the receptor conversion rate βr). The neural network takes as input two vectors,
, representing expression of all receptors in cell j, and
, representing the average expression of all ligands across the set of all N cells. These two vectors are then fed into the neural network (F) to predict the expression
of all genes in cell j (regression task). The neural network structure comprises the following components: two input layers, two hidden layers, and one output layer. The size of the two input layers corresponds to the number of ligands and receptors, respectively. The first hidden layer is dedicated to representing specific LR pairs. Each node in this layer receives two connections, one corresponding to a ligand and the other to a receptor, thus encapsulating the unique LR pair interactions. The second hidden layer represents transcription factors involved in downstream signaling pathways. We leverage scSeqComm (Baruzzo et al. 2022), a database that incorporates insights from KEGG (Kanehisa and Goto 2000) and REACTOME (Jassal et al. 2020) pathways to delineate the association between a receptor and its downstream transcription factors, to inform the connections between this layer and the previous one. The output layer is responsible for outputting downstream gene expression predictions. To establish connections between known TFs and their downstream gene targets, we utilize prior knowledge from databases such as TRRUST v2 (Han et al. 2018), HTRIdb (Bovolenta et al. 2012), and RegNetwork. If specific information regarding the target genes of certain TFs is unavailable, we utilize dense connections. We use a stochastic gradient descent optimizer to train the network F that minimizes the mean squared error loss L (between the true observed expression for the set of all cells X and the predicted expression):
(2)
With the trained neural network regressor, we then perform a sensitivity analysis to score each input receptor based on its impact on the downstream gene expression. Specifically, we conduct a perturbation-based feature calculation (Mi et al. 2021), by permuting feature values for each receptor, independent of the output. Next, we adjust the expression of the permutated receptor to reflect a ∼50% reduction, scaling it by the median expression of all receptors. This represents the perturbed input receptor vector (
) simulating a ∼50% repression of receptor r in expression. Although this 50% repression serves as the default setting, we also offer the option to simulate a 100% repression, equivalent to a complete knockout, catering to users who prefer a full knockout scenario. Then, we quantify the importance of the receptor r (Ir) based on the loss increase caused by the perturbation:
(3)
Employing a similar approach as the receiving rate, we normalize the receptor importance into range (0,1) via min-max normalization. Here, we set a minimum value for the conversion rate to be 0.4 based on previous studies that measure the effect of cellular interaction on gene expression (Arnol et al. 2019). We confirmed that this approach does not significantly alter the resulting gene expression of the cells of the training grid (Supplemental Figs. 28, 29).
We included various validations of our deep learning module to estimate the conversion rate as well its capabilities to perform in silico perturbation (see Supplemental Results; Supplemental Figs. 30–33). We also include ablation studies that show the positive impact of using biologically informed neural network using a gene network and the conversion rate as a whole, in line with other studies (Supplemental Fig. 34; Lin et al. 2017). Please refer to the Supplemental Results for further details.
-
Cell–cell interaction score between cells: The interaction score between a ligand–receptor pair (l, r) in a sender-receiver pair (i, j) is calculated based on the rates calculated above (γl, αr, βr):
(4)
The cell receiving score for a given cell i is calculated as the summation over all ligand–receptor pairs (l,r) across all sending cells:
(5)
The combined interaction score between a sender-receiver cell pair (i, j) is calculated as the summation over all ligand–receptor pairs (l,r):
(6)
The overall interaction score IS(C1, C2) between two clusters C1 (sender) and C2 (receiver) for a ligand–receptor pair (l,r) is quantified as the average interaction score of all possible sender-receptor cell agent pairs between the sender C1 and receiver C2:
(7)
We also introduce a method to optionally integrate the utilization of pseudotime ordering of cells into the calculation of CCI scores. Please refer to the Supplemental Methods for further details.
Unlike typical ABMs, CellAgentChat integrates various statistical methods to estimate parameters governing agent behaviors. However, it stands out as a robust ABM framework rather than a workflow. Whereas workflows consist of discrete, sequential steps operating independently or semi-independently, frameworks provide an integrated structure where components function cohesively as a unified system. CellAgentChat adheres to core ABM principles (Bonabeau 2002; Masad and Kazil 2015; Yu and Bagheri 2020), modeling cells as agents with unique states (e.g., gene expression profiles and spatial coordinates) that dynamically evolve through interactions with their environment and neighboring cells. The statistical and neural network components in CellAgentChat are not standalone modules, as in a workflow, but integral elements that enhance the framework's ability to simulate realistic interaction rules from data.
Finding significant ligand–receptor pairs between cell populations
After computing interaction scores (IS) for LR pairs between cell clusters, we assess their statistical significance using
a random permutation test. Specifically, we randomly permute cell group labels and recalculate LR interaction scores, repeating
this until each LR pair accumulates approximately 10,000 random scores. We then estimate the parameters (alpha and scale)
of a gamma distribution from these scores. The significance (P-value) of each observed IS is computed based on its position in the right tail of this distribution, followed by FDR (Benjamini–Hochberg)
correction. Interaction scores with adjusted P-values below 0.05 are deemed significant (Significant Interaction Scores, SIS). Using the SIS, we can calculate the cell
type pair scores between two clusters, C1 (sender) and C2 (receiver):
(8)
Because of the computational expense of calculating the distances between all cells in the data set, the permutation test
can be quite costly when considering spatial information. To address this, we have devised a method to match the background
gamma distribution of non-spatial calculations with that of spatial ones. This is achieved by dividing the scale parameter
by the average distance between all cells. Additionally, we halve both the alpha and the modified scale values. Although this
does not equate to perfect scaling, we present very similar results using this method significantly more efficiently (Supplemental Fig. 35B,C).
In silico receptor perturbation
After identifying the top interactions between cell types, we conducted in silico perturbation using our trained deep learning model F. To illustrate the process of in silico receptor blocking, we centered our attention on the receptors involved in the top 25 LR interactions, spanning all cell type pairs. We then subjected these receptors to a perturbation analysis.
We permutated feature values for each receptor, independent of the output and simulated receptor blocking and scaled the expression
of the permutated receptor by the median expression of all receptors. This manipulation of the input receptor vector, , allowed us to gauge the ensuing influence on the gene expression downstream. After each receptor's perturbation, we singled
out the most impacted genes, on average, across all cells. Specifically, we looked at the top 50 most impacted genes to assess
the effects of receptor blocking at the transcriptome level. To further demonstrate the potency of our model in executing
in silico perturbation, we utilized the DisGeNET (Piñero et al. 2015) database as the potential validation. We verified whether the perturbed genes are associated with the specific diseases
that piqued our interest—breast cancer/PDAC (i.e., diseases associated with the data set). To quantify the strength of the
association between each receptor and the disease, we applied a binomial test to compute the statistical significance of the
overlap between the perturbed genes and the genes associated with each disease (breast cancer, PDAC). We further refined our
results by adjusting the P-values over all disease gene sets, utilizing the FDR-Benjamini–Hochberg method, ensuring the robustness of our findings.
Additionally, we confirmed the effectiveness of these receptor-inhibitor therapeutics using Kaplan-Meier survival plots (Győrffy 2021).
To validate CellAgentChat's in silico receptor blocking functionality, we compared its performance with scGen (Lotfollahi et al. 2019). We trained scGen to reconstruct cell gene expression profiles and reduced each receptor's expression by 50% during inference to replicate CellAgentChat's strategy. The model predicted the resulting downstream gene expression changes, and we analyzed these changes to identify the most affected genes. We then cross-referenced these genes with DisGeNET to evaluate their relevance to the disease of interest, using the same evaluation pipeline as CellAgentChat.
Agent-based dynamic simulations and animation of CCIs
CellAgentChat offers dynamic simulation capabilities for analyzing the impact of CCIs on cell states over short durations. Users can either execute the model once, similar to existing methods, inferring CCIs at a single moment, or run it across multiple short time steps to simulate temporal dynamics. At each step, CCIs are computed, and a trained neural network predicts gene expression levels for the next step. Although predicted gene expression values cannot directly replace actual data owing to uncertainties, we trust their directional changes—whether genes increase or decrease (Ding et al. 2018). Genes with statistically significant changes are then adjusted slightly (±0.01) within the ABM, simulating minor expression variations. This iterative process reflects the dynamic modulation of gene expression by CCIs.
Uniquely, CellAgentChat infers interactions between individual cells, not just cell types. Its real-time animation platform visualizes cell response scores, where size indicates CRS strength—larger cells represent stronger responses (Supplemental Fig. 35A). Users can customize hyperparameters, such as the number of steps (default = 1), δ (default = 1), and τ (default = 2), to perform effective in silico permutations. See Supplemental Methods for more information to effectively set up hyperparameters to facilitate simulations.
In addition, CellAgentChat offers hyperparameters for controlling the LR universe. Users can block a specific LR pair or a receptor to observe downstream gene expression changes. The animation platform also allows for visualizing specific cell types, with modifications only applied to those selected. Furthermore, users have the option to visualize specific receptors or LR pairs.
Further details on standard data set preprocessing and methods related to benchmarking and functional evaluation of inferred ligand–receptor pairs can be found in the Supplemental Methods.
Data sets
CellTalkDB is available at GitHub (https://github.com/ZJUFanLab/CellTalkDB). The SCC and PDAC data sets analyzed in this study are available from the NCBI Gene Expression Omnibus (GEO; https://www.ncbi.nlm.nih.gov/geo/) repository under the following accession numbers: GSE144240 (GSM4565826) and GSE111672 (GSM3036911), respectively. The data for the SeqFISH+ mouse somatosensory cortex data set can be found at GitHub (https://github.com/CaiGroup/seqFISH-PLUS). The Xenium human breast cancer data set (sample 2) can be found at https://www.10xgenomics.com/products/xenium-in-situ/preview-dataset-human-breast. The Stereo-seq developmental mouse hippocampus data set can be found at STOmicsDB (Xu et al. 2024) https://db.cngb.org/stomics/mosta/. We used the h5ad file titled “Dorsal_midbrain_cell_bin.h5ad”. The PDAC data set with tumor and healthy regions can be found at Zenodo (https://zenodo.org/records/10712047).
Software availability
The CellAgentChat software package and source code are available at GitHub (https://github.com/mcgilldinglab/CellAgentChat) and as Supplemental Code.
Competing interest statement
The authors declare no competing interests.
Acknowledgments
This work is supported by grants from the Canadian Institutes of Health Research (CIHR) (PJT-180505 to J.D.); the Fonds de recherche du Québec - Santé (FRQS) (295298 and 295299 to J.D.); the Natural Sciences and Engineering Research Council of Canada (NSERC) (RGPIN2022-04399 to J.D.; RGPIN-2016-05174 to Y.L.); and the Meakins-Christie Chair in Respiratory Research (to J.D.). This publication is part of the Human Cell Atlas (HCA) publication bundle (HCA-10).
Author contributions: J.D. conceived the study. J.D. and Y.L. supervised the study. V.R. built the framework. V.R. conducted the experiments, with assistance from Y.Z. on one experiment. V.R., J.D., and Y.L. wrote and revised the manuscript. All authors read and approved the final manuscript.
Footnotes
-
[Supplemental material is available for this article.]
-
Article published online before print. Article, supplemental material, and publication date are at https://www.genome.org/cgi/doi/10.1101/gr.279771.124.
-
Freely available online through the Genome Research Open Access option.
- Received July 9, 2024.
- Accepted April 23, 2025.
This article, published in Genome Research, is available under a Creative Commons License (Attribution-NonCommercial 4.0 International), as described at http://creativecommons.org/licenses/by-nc/4.0/.