PearlST comprises three integral components: gene expression augmentation, low-dimensional representation learning, and downstream analysis (Fig. 1). First, PearlST enhances gene expression by incorporating information from histology images, spatial location, and gene expression profiles of spatially nearest neighbors for each spot. Notably, PearlST leverages ResNet-50 [29] within the SimCLR framework [30] to learn visual features of each spot image, utilizing histology images collected from diverse tissues across species. To ensure spatial proximity consistency across histology images, spatial locations, and gene expression, the 3 datasets are integrated to identify spatial nearest neighbors for each spot. It is important to emphasize that the feature extraction module is exclusively employed for sequencing platforms that generate histology images, such as 10X Visium. For other ST data, spatial coordinates and gene expression are used to compute the nearest neighbors of the spots. Subsequently, gene expression data undergo denoising and enhancement through a PDE-based diffusion model. Secondly, PearlST learns latent low-dimensional embeddings that encompass both spatial information and augmented gene expressions by employing a Wasserstein adversarial regularized graph autoencoder. The introduction of a Wasserstein regularizer minimizes the 1-Wasserstein distance [31] between distributions. Thirdly, the obtained latent embeddings are utilized for data visualization through Uniform Manifold Approximation and Projection (UMAP) [32]. Based on the latent embeddings, spatial domains are segmented using common clustering algorithms, defaulting to K-means or, alternatively, Louvain clustering. Furthermore, a pSM [27] of the tissue is computed by employing the diffusion pseudotime (DPT) method [33] based on the latent embeddings (Fig. 1).

To assess the efficacy of PearlST in identifying distinct layer-specific expression profiles, we initially applied it to 12 human DLPFC samples. These samples were previously annotated into 6 cortical layers and white matter (WM) by Maynard et al. [34], considering cytoarchitecture and selected gene markers. Using the manual annotation results as the ground truth, we employed the adjusted Rand index (ARI) to quantify the similarity between the inferred domains and the expert annotations across all 12 sections. In this evaluation, we compared the performance of PearlST with 8 recently proposed spatial clustering approaches (BayesSpace [19], DeepST [25], SpaGCN [26], stLearn [24], SEDR [22], SpaceFlow [27], STAGATE [23], and SpecMix [21]), as well as 2 non-spatial clustering methods implemented by SCANPY [35] and Seurat [36]. We ran each method according to the official technical documentation and used the default parameters for each method (Note S1). The comparison results reveal that PearlST achieves a significantly higher median ARI score than the other methods across the 12 samples (Fig. 2A). Following that, STAGATE, DeepST, and BayesSpace have higher median ARI scores on these datasets. We conducted a one-sided Wilcoxon rank-sum test comparing the results of PearlST with these 3 methods. The statistical test results indicate that PearlST performs significantly better in clustering on the DLPFC dataset compared to STAGATE (P = 6.1 × 10⁻³), DeepST (P = 7 × 10⁻⁴), and BayesSpace (P = 2 × 10⁻⁴). Notably, SEDR and SpecMix exhibit comparable clustering performance.

Subsequently, we conducted a more detailed analysis on section 151671 (Fig. 2B to D). We visualized the domain segmentation results obtained by each method and compared them to expert annotations (Fig. 2B). Notably, all methods exhibited challenges in capturing the subtle structure of Layer 4, potentially attributed to the spatial resolution limitations inherent in this specific ST data. Despite this, both PearlST and BayesSpace successfully captured the remaining structures (Layer 3, Layer 5, Layer 6, and WM) as observed in the manual annotations. Furthermore, it was observed that the expression profiles of specific layers segmented by BayesSpace had more pronounced jagged edges, while the results from PearlST appeared relatively smoother. Specifically, for the subtle Layer 4 structure, BayesSpace annotated it as a large area, whereas PearlST provided a smaller and more refined annotation result. Other methods also seemed to dissect section 151671 into 5 layers, but there was more noise at the edges of layer-specific expression profiles, particularly for SCANPY and Seurat. Similar comparison results were observed across the remaining 11 samples (Figs. S1 to S3).

Additionally, we illustrate that PearlST excels in revealing distances between spatial domains and depicting spatial trajectories in UMAP plots. To achieve this, we utilized SCANPY to reduce the low-dimensional representation of the output from various methods to 2 dimensions for UMAP visualization and PAGA (partition-based graph abstraction) trajectory inference. Notably, the UMAP plots generated by embeddings from PearlST and STAGATE exhibited well-organized cortical layers and consistent spatial trajectories, progressing from WM to Layer 3 (Fig. 2C and Fig. S4). Among the other methods, only DeepST and SpecMix results demonstrated somewhat similar consistency to PearlST. In the PAGA graphs (Fig. S5), the low-dimensional representations of PearlST displayed an approximately linear developmental trajectory from Layer 3 to Layer 6, with stronger similarity observed between adjacent layers compared to non-adjacent layers. Although the low-dimensional representation of DeepST also showcased a linear trajectory, it did not correctly depict the progression from Layer 3 to Layer 6.

In the final analysis, we compared the pSM inferred by different methods, specifically focusing on spatially aware methodologies. STAGATE, SpaceFlow, SEDR, and PearlST all presented layer-patterned pSMs with a distinct and smooth color gradient (Fig. 2D). Notably, these maps suggested a pseudo-spatiotemporal ordering from WM to Layer 3, aligning with the developmental order of the cortex from the inside out and reflecting the hierarchical spatial organization of the tissue. However, it is crucial to highlight that only the results from PearlST exhibited a clear and unique developmental starting point, along with a smooth progression between neighboring layers (Fig. 2D and Fig. S6).

Subsequently, we illustrate the efficacy of the PDE-based diffusion model in noise reduction for ST data, emphasizing its ability to characterize spatial expression patterns. PearlST was applied to the 12 DLPFC samples utilized in the previous section, and ablation experiments were designed to investigate the impact of data enhancement on downstream analyses. The findings revealed that enhancing gene expression data significantly improved the accuracy of spatial domain identification (Fig. S7). For further validation, we selected 5 of these samples to demonstrate that the augmented data facilitated other spatiotemporal structure inference tasks as well (Fig. S8). In the low-dimensional visualization results, the enhanced data showcased a clear hierarchical structure in samples 151507, 151670, 151672, and 151675, while the results from the raw data were accompanied by a substantial amount of noise, and the original hierarchical structures were not even discernible in samples 151507 and 151672 (Fig. S8A). The enhanced gene expression data notably contributed to the generation of more consistent PAGA trajectory inferences, particularly evident in samples 151670 and 151672 (Fig. S8B). On these few samples, neither the raw nor the augmented data produced pSMs that exhibited a high degree of agreement with the true developmental order; however, the latter results demonstrated a clearer and more progressive order (Fig. S8C).

Furthermore, we applied PearlST to an image-based ST dataset generated by STARmap [14] with single-cell resolution. The STARmap dataset comprises 1,020 genes expressed on 1,207 cells and has been annotated with 7 structural layers (Fig. 3A) and 16 cell types (Fig. 3B). Similar to previous analyses, we initially clustered the ST data using PearlST, considering expert annotations as the gold standard, and compared the results with 4 other methods. The number of clusters for all compared approaches was set to 6. PearlST demonstrated the most accurate domain segmentation with an ARI score of 0.67. SpaceFlow followed with an ARI of 0.60 and STAGATE was 0.54, while the remaining 3 methods yielded clustering results with ARI scores below 0.4 (Fig. 3C). These results are further manifested in the distribution of specific regions, revealing more overlap between regions produced by SCANPY and SEDR. In contrast, SpaceFlow, STAGATE, and PearlST exhibited significantly less overlap. While SpaceFlow showed slight overlap on CC and HPC and larger overlap on HPC and L6, wrongly identifying a portion of L1 as L2/3, PearlST only annotated a small portion of L1 to L2/3 and very few to L5 (Fig. 3C). The outputted low-dimensional representations from these methods were used to compute the pSM. SCANPY exhibited the least satisfactory result, lacking a clear differentiation trajectory (Fig. 3G). Both STAGATE and SEDR results did not show a gradient process, with dark green dominating the former and light green dominating the latter (Fig. 3E and F). SpaceFlow demonstrated a satisfactory progression from the left to the center; however, it presented an anomalous jump from the very first developmental stage on the right to the final form in the center (Fig. 3D). In contrast, the results from PearlST were the most satisfying, displaying a gradual color change from darker to lighter from left to right, effectively unveiling the developmental trajectory of the mouse visual cortex (Fig. 3H).

We demonstrate here that PearlST excels in delineating the laminar structure of data acquired with Stereo-seq [37]. Stereo-seq is an emerging spatial omics technology that provides subcellular spatial resolution through DNA nanopore pattern array chips. For validation, we utilized a coronal MOB tissue dataset, a widely used model tissue with laminar organization. Fu et al. [22] annotated the laminar organization in the DAPI-stained image, containing the rostral migratory stream (RMS), granule cell layer (GCL), internal plexiform layer (IPL), mitral cell layer (MCL), external plexiform layer (EPL), olfactory nerve layer (ONL), and glomerular layer (GL) (Fig. 4A). SpaceFlow and PearlST partitioned this dataset into 7 layers, and STAGATE annotated the dataset into 8 layers according to original tutorial (Fig. 4B). Both PearlST and STAGATE provided well-defined layer structures and effectively demarcated the GL and MCL regions, while SpaceFlow results appeared to have extensive overlap. Subsequently, we applied marker genes from each anatomical region to validate PearlST results (Fig. 4C and D). PearlST demonstrated good correspondence between clusters and known marker genes. Some marker genes, such as Mbp and Pcp4, showed high expression levels overlapping with neighboring regions, which is expected due to shared cell types in different internal structures of organs. Overall, PearlST accurately identified relevant anatomical regions. Finally, we computed the pSM using the low-dimensional representations of each method. The pSM values (green) were lowest in the EPL and increase as they move toward the interior. This sequence of pSM coincides with the developmental order of these layers, starting in the central EPL and developing outward toward the sides, leading to the MCL and the GL, ONL, and GCL developing last [38]. Values are highest in the interior, the peaks being in the GCL and RMS. This suggests that spatial flow calculations of pSM are not only clearer than non-spatial pseudotime estimation but also more accurately reflect the temporal and developmental relationships between cells. Overall, the pSM calculated by PearlST (Fig. 4G) is in high agreement with the real developmental order in these layers. In contrast, the pSM results from both STAGATE (Fig. 4E) and SpaceFlow (Fig. 4F) indicate a developmental direction from the outermost IPL layer to the middle-most MCL layer, which is inconsistent with the true developmental order.

In this example, we demonstrate that PearlST excels in revealing finer-grained tissue structures using Stereo-seq data of a mouse embryo at E12.5, containing 55,295 bins and 27,330 genes [37]. Chen et al. [37] annotated the mouse embryo E12.5 transcriptional profile into 19 domains (Fig. 4H). Initially, we performed spatial clustering of this dataset using PearlST, setting the number of clusters to 19. The clustering results by PearlST matched the annotated regions well (Fig. 4I). PearlST effectively captured most of the fine-grained structures in the embryo, demonstrating high consistency with the original annotations. Major regions such as the cartilage primordium, brain, heart, dorsal root ganglion, liver, and muscle regions were accurately identified (Fig. 4I). However, there were 3 regions (blood vessel, lung, and sympathetic nerve) that PearlST did not annotate accurately. Further visualization and comparison of annotations by PearlST and the ground truth for each major region demonstrated higher consistency between them (Fig. 4J and K). Importantly, PearlST annotations exhibited high concordance with known marker genes of major organs (Fig. 4L). For instance, the cartilage primordium region is marked by Meox1, Pax1, and Pax9; the dorsal root ganglion region is marked by Prdm12 and Ntrk1; epidermis is marked by Krt5, Krt14, and Krt15; heart is marked by Nppa, Myl7, Myl2, and Myh7; liver is marked by Afp, Alb, and Fgb (Fig. 4L and Figs. S10 and S11). Additionally, there are regions marked by individual genes, such as the brain region marked by Tubb2b and the meninges region marked by Map1b.

We further employed Slide-seqV2 ST data with a spatial resolution of 10 μm to evaluate the performance of PearlST. Initially, we applied PearlST to a Slide-seqV2 dataset from the mouse hippocampus [7], annotated by the Allen Reference Atlas [39] (Fig. 5A). As anticipated, PearlST effectively characterized the tissue structure, revealing distinct spatial domains. In comparison, the clusters identified by SpaceFlow lacked significant spatial separation (Fig. 5B). Notably, PearlST and STAGATE depicted a clear “rope-like” structure and an “arrowhead-like” structure in the hippocampal region, accurately identifying its 5 spatial domains. This result is consistent with the annotation of hippocampal structures in the Allen Reference Atlas (Fig. 5B). All 3 methods used for comparison produced spatially consistent clustering, capturing major anatomical regions such as the dentate gyrus (DG) and the pyramidal layers within Ammon's horn, further divided into CA1, CA2, and CA3 regions. In particular, PearlST and STAGATE exhibited sharper boundaries in delineating the CA3 and DG regions, aligning more closely with anatomical annotations (Fig. 5A and C) and marker gene expression (Fig. 5D). In contrast, SpaceFlow annotated CA1, CA3, and DG as the same region. Additionally, both PearlST and SpaceFlow annotated the LH and V3 regions. To validate the clustering performance of PearlST, we examined the marker gene expression of the regions (Fig. 5D), revealing good agreement between most regions and their corresponding marker genes. This further supports the robustness and accuracy of PearlST in ST data analysis.

Furthermore, we applied PearlST to sections of the MOB sequenced by Slide-seqV2 [7] and observed that the spatial domains identified by PearlST were consistent with the annotations in the Allen Reference Atlas (Fig. 5E). As previously mentioned, Fu et al. [22] annotated the laminar organization of the coronal MOB in the DAPI-stained image, including the RMS, GCL, IPL, MCL, EPL, ONL, and GL. In this case, PearlST annotated the dataset into 10 regions, as in a previous study [23]. Specifically, PearlST identified 2 spatial domains corresponding to the accessory olfactory bulb (AOB) and the granular layer of the accessory olfactory bulb (AOBgr), respectively (Fig. 5F).

The spatial domains uncovered by PearlST are strongly supported by known gene markers (Fig. 5G and H). For instance, the gene Nrgn exhibited high expression in the GCL_1 domain, while the gene Gap43 showed strong expression in the segmented region AOB, consistent with immunohistochemical experimental results [40]. The granular cell marker Atp2b4 [41] displayed strong expression in the identified AOBgr domain, and the mitral cell marker Gabra1 [42] exhibited dominant expression in the narrow MCL structure, also identified by PearlST. Notably, PearlST identified a spatial subgroup called GCL_1 that is predominantly expressed by Nrgn. Considering that Nrgn is a well-documented risk gene for schizophrenia [43], this suggests a potential association between this structural domain and cognitive function. In summary, these results demonstrate the capability of PearlST to identify spatial domains in Slide-seqV2 sequencing data with much higher spatial resolution, enabling the exploration of tissue structures in finer detail.

To investigate the spatial heterogeneity of the cancer microenvironment, we applied PearlST to a 10X Visium ST dataset of human breast cancer. Initially, we identified spatial domains for the human breast cancer dataset and compared the results with manual annotation based on pathological labeling from H&E staining by SEDR [22]. Notably, unlike previous datasets on DLPFC, mouse hippocampus, etc., which have clear and defined morphological boundaries, tumor tissues are highly heterogeneous and contain a complex microenvironment [22,44]. Based on pathological features, the histological image was segmented into 20 regions, consisting of 4 main morphotypes: ductal carcinoma in situ/lobular carcinoma in situ (DCIS/LCIS), healthy tissue (Healthy), invasive ductal carcinoma (IDC), and tumor surrounding regions with low features of malignancy (Tumor edge) (Fig. 6A). We assessed the similarity between clustered results and ground truth using the ARI score, comparing a range of popular spatial clustering methods. For Seurat and SCANPY, many of the computed clusters appeared fragmented and discontinuous (Fig. 6D), while STAGATE, SpaGCN, SpaceFlow, DeepST, and PearlST produced fewer disjointed clusters (Fig. 6B to D). PearlST exhibited the best results in terms of region segmentation accuracy and edge noise, particularly for the identification of Tumor_edge_2, DCIS_LCIS_4, IDC_2, IDC_4, IDC_5, IDC_6, and Healthy_1. These results were also reflected in the ARI scores, with PearlST achieving the highest score at 0.65 (Fig. 6B), while DeepST and stLearn scored under 0.55, and SEDR, SpaceFlow, and SCANPY scored around 0.49. The remaining 3 methods had lower ARI scores (Fig. 6C and D).

Continuing our analysis, we extended PearlST to study cell–cell communication between regions. Initially, we selected 9 regions that were closer to the manual annotations to explore the strength of their interactions (Fig. 6E). On a macro level, stronger interactions were observed between regions DCIS/LCIS_4 and IDC_2, as well as between DCIS/LCIS_4 and IDC_4. Tumor_edge_2 received a significant number of ligand-signaling molecules from IDC_4 and IDC_8. Notably, weaker interaction strengths were observed between regions of the same type, such as IDC, Tumor edge, and DCIS/LCIS (Fig. 6E). Zooming in on a micro level, we studied a multilayer signaling network between any 2 regions. Taking the example of constructing a multilayer communication network between 2 regions IDC_4 and Tumor_edge_1, we employed the law of mass reciprocity to calculate the ligand–receptor (LR) scores between these regions. We then used random forest regression to establish the relationship between LR scores and low-dimensional embeddings of spots (Fig. S12A and B). To further explore the genes contributing significantly to the low-dimensional embeddings, we computed feature importance using stochastic perturbations. Genes identified as significantly contributing genes were considered as target genes in the multilayer network. The role of LR pairs in regulating target genes was reflected by their previous effects on low-dimensional embeddings. The analysis revealed that certain LR pairs, including TNF-TNFRSF21, RBP4-STRA6, PDGFA-PDGFRB, TNC-CNTN1, and ALKAL2-ALK, exhibited stronger interactions (Fig. 6F). In particular, ALKAL2-ALK exerted a dominant regulatory effect on the generation of embedding_9, while PDGFA-PDGFRB regulated embedding_5 and embedding_12 simultaneously (Fig. 6F and Fig. S12C to G).

Lastly, we performed functional enrichment analysis of target genes using clusterProfiler [45], encompassing KEGG (Kyoto Encyclopedia of Genes and Genomes) enrichment (Fig. S13A) and GO (Gene Ontology) enrichment (Fig. S13B). The analysis revealed that the target genes are significantly enriched in numerous signaling pathways associated with cancer development and progression. Notable pathways include the PI3K-AKT signaling pathway [46], p53 signaling pathway [47], immunoglobulin-mediated immune response [48], and B cell receptor (BCR) signaling pathway [49] (Fig. S13). Notably, the abovementioned ALKAL2-ALK signaling [50–53] and PDGFA-PDGFRB signaling [54,55] have been reportedly important in breast cancer growth and metastasis. Moreover, the PI3K-AKT pathway has been validated as the downstream signaling pathways of ALKAL2-ALK signaling and PDGFA-PDGFRB signaling in human cancers [53,56,57]. These findings suggest important functional interactions between the 2 regions IDC_4 and Tumor_edge_1 through LR interactions in the tumor microenvironment. More specifically, from the multilayer and functional regulation network of the learned embedding, one can infer which target genes are important genes contributing to this embedding and how they are regulated by other domains via intercellular and intracellular signaling pathways.

The preprocessing of gene expression matrix involves the following steps. Initially, gene filtering is executed to retain only those genes expressed in a minimum of 4 cells or more. Subsequently, normalization occurs, wherein the expression of each gene is divided by the total expression in the respective cell. This normalization ensures that each cell attains an identical total count after the process. The normalized expression is then multiplied by a scale factor (default 10,000) and subjected to log-transformation with a pseudo-count of one. Lastly, the top 2,000 highly variable genes are selected as the inputs for PearlST. The above procedures were implemented using the SCANPY package.

Spatial gene expression technology offers transcriptome-wide gene expression profiles, accompanied by spatial position information and/or tissue morphology. In the case of PearlST, these additional tissue data are leveraged to enhance gene expression in adjacent spots. PearlST employs ResNet-50 [29] trained through the SimCLR framework [30] to learn visual features from each spot image. The training dataset comprises histology images collected from various tissues of different species. PearlST evaluates gene expression similarity, morphological similarity, and spatial relationships between each pair of spots. The gene expression weight between spot i and spot j (ES_ij) was computed using the cosine similarity as follows:where G_i and {\overline{G}}_i denote the expression profiler and the mean expression of spot i, respectively. Following this, the morphological similarity between spot i and the adjacent spot j (MS_ij) was computed as follows:where M_i refers to the visual features of spot image i, with each spot image transformed into 2,048 dimensional latent variables. Subsequently, the spatial distance between the considered cell and all other cells is calculated. The distances between the top 30 (optional) adjacent spots are then ordered to determine the radius γ (mean add variance). For a given spot i, a spot S_j is considered to be a neighbor, then PS_ij = 1 if and only if the Euclid distance between the 2 spots is less than γ; otherwise, PS_ij = 0. Finally, under the combined consideration of spatial gene expression, histological similarity, and spatial adjacency, the overall pseudo similarity (OS) between spots i and j can be expressed as:We denote by NN_i the set consisting of the top (default by 4) nearest neighbors of spot i according to Eq. 3. In this way, we integrate multiple pieces of information to more precisely calculate the spatial nearest neighbors of each cell or spot.

Given the inherent noise in spatial transcriptome data, we utilized a PDE-based diffusion model to both denoise and enhance gene expression. This enhancement is achieved by considering the spatial nearest neighbors of each spot, as obtained in the previous section. Here, the transcriptional profile of log-transformed gene expression for each spot is represented by u(x,  y) = u(x,  y,  t)|_t=0 in a bounded region Ω. Then, the following PDE model is employed to describe the diffusion process of u(x,  y),where c(|∇u|) is referred to as the diffusion function, which is a function of the gradient mode of the image. In other words, the smoothing operator of this diffusion model varies with the changes in the magnitude of the gradient. For this purpose, we adopt a diffusion function proposed by Perona and Malik [59]:where K is a constant, taken by default to be 0.2. Discretizing Eq. 4 gives its iterative difference form as:In the above equation, \nabla S{u}_{i,j}^n=u_{i+1,j}^n-u_{i,j}^n, \nabla E{u}_{i,j}^n=u_{i,j+1}^n-u_{i,j}^n, \nabla N{u}_{i-1,j}^n=u_{i-1,j}^n-u_{i,j}^n, \nabla W{u}_{i,j}^n=u_{i,j-1}^n-u_{i,j}^n, and λ is a small positive parameter, which, in this article, defaults to 0.1. See Note S3 for a more detailed derivation. Prior to data augmentation, the transcriptional profile of each spot is reshaped to a 50 × 40 gray scale plot. Gene expression denoising is achieved by iteratively updating Eq. 6, where u_{i,j}^n denotes the pixel value at grid (i, j) and n represents iteration times. Next, gene expression augmentation is performed by substituting the gene expression of the corresponding spatial nearest neighbors of each spot obtained by Eq. 3 into Eq. 6 for iteration. This means that \nabla S{u}_{i,j}^n, \nabla E{u}_{i,j}^n, \nabla N{u}_{i-1,j}^n, and \nabla W{u}_{i,j}^n are computed using the expressions of the nearest neighbors belonging to NN_i. The number of iterations is set to 4 in this study, which reconciles both prediction performance and computational efficiency.

The domain segmentations are acquired by applying K-means to the low-dimensional embeddings from PearlST. The pSM is computed by employing the DPT algorithm [33], utilizing the low-dimensional representations output from PearlST. DPT employs diffusion-like random walks to estimate the ordering and transitions between cells. By utilizing the embeddings from PearlST, which encode both spatial location and gene expression profiles of cells, DPT generates a spatiotemporal order that is consistent in both space and pseudotime. The root cell of the pSM can be specified with a priori knowledge; otherwise, in our strategy, the cell in the embedding space that has the largest sum distance from the other cells is assigned as the root cell.

Building on our prior work stMLnet [62], we quantify spatially dependent LR signaling activity and infer downstream functions by associating LR pairs to the latent features learned by PearlST and target genes. Specifically, the spatial distance-dependent LR signaling scores can be quantified based on a diffusion model and the law of mass action. After obtaining LR scores, we employed Random Forest regression to learn the regulatory relationship between LR pairs and the low-dimensional representation of receiver cells. Additionally, we computed the feature importance of the graph autoencoder by random perturbations identifying those genes that contribute more to the low-dimensional representations and considering them as target genes. This process allowed us to construct a multilayered signaling network (ligand–receptor–embeddings–target genes) for dissecting cellular communications between different spatial domains. Furthermore, enrichment analysis (e.g., KEGG or GO enrichment) of target genes can be performed using tools like clusterProfiler [45] to elucidate the biological functions of the LR-mediated cellular interactions. For a more detailed description of the inference of functional and multilayer cellular communication networks, please refer to Note S4.