Screening cell–cell communication in spatial transcriptomics via collective optimal transport

Spatial transcriptomic technologies and spatially annotated single-cell RNA sequencing datasets provide unprecedented opportunities to dissect cell–cell communication (CCC). However, incorporation of the spatial information and complex biochemical processes required in the reconstruction of CCC remains a major challenge. Here, we present COMMOT (COMMunication analysis by Optimal Transport) to infer CCC in spatial transcriptomics, which accounts for the competition between different ligand and receptor species as well as spatial distances between cells. A collective optimal transport method is developed to handle complex molecular interactions and spatial constraints. Furthermore, we introduce downstream analysis tools to infer spatial signaling directionality and genes regulated by signaling using machine learning models. We apply COMMOT to simulation data and eight spatial datasets acquired with five different technologies to show its effectiveness and robustness in identifying spatial CCC in data with varying spatial resolutions and gene coverages. Finally, COMMOT identifies new CCCs during skin morphogenesis in a case study of human epidermal development.


Supplementary Figure 9
Inferred sent or received signal exceeding the amount of ligand or receptors Each dot represents a ligand-receptor pair. The commonly used KL divergence-based unbalanced optimal transport is used. All box plots show Q1 (25% quantile), Q2 (median), Q3 (75% quantile) as the box and the upper whisker represents Q3+1.5IQR. All data points are 0% for every COT experiments demonstrating no exceeded mass. For the Visium human breast cancer and Visium mouse brain data, the boxplots contain 1026 and 960 ligand-receptor pairs, respectively.

Supplementary Figure 10
Comparison of whether to normalize expression to probability distributions.
a An example with 8 cells where 1 cell expresses the ligand gene and the other 7 cells express the receptor gene. b The inferred cell-cell communication (CCC) with or without normalizing the expression levels to probability distributions. c An example ligand-receptor pair in a Visium human breast cancer data. d The inferred spatial distribution of received signals through the ligandreceptor pair FGF19-FGFR1 with or without normalization to probability distribution and with or without a spatial distance cutoff.

Supplementary Figure 11
Comparison of whether to incorporate spatial distance information.
a An example with 8 cells where 1 cell expresses the ligand gene and the other 7 cells express the receptor gene. b The inferred cell-cell communication (CCC) with or without normalizing the expression levels to probability distributions. Both results are from the spatial expression levels without normalizing to probability distributions. c An example ligand-receptor pair in a Visium human breast cancer data. d The inferred spatial distribution of received signals through the ligandreceptor pair FGF19-FGFR1 with or without incorporating spatial information and if spatial information is used, with or without a spatial distance cutoff.

Supplementary Figure 12
Leave-one-out cross validation of scRNA-seq and spatial data integration for the human epidermis data For each gene in the spatial data, the integration of scRNA-seq and spatial data is performed using the rest of the genes and the predicted spatial expression of this gene is compared to the groundtruth.

Supplementary Figure 15
JAK-STAT signaling in human skin a Cell-cell communication directions. Sender: the direction shows to where the signal sending cells are sending the signal and color depicts the amount of signal sent by each cell. Receiver: the direction shows from where the signaling receiving cells are receiving the signal and color depicts the amount of signal received by each cell. b The cell-cell communication of all ligand-receptor pairs summarized to clusters (clustering identical to the recent NC paper), a thicker arrow means stronger signaling activity. c Dotplot showing the detailed signaling strength among the clusters through the individual ligand-receptor pairs. d The identified differentially expressed genes due to cell-cell communication activity. This is analogous to pseudotime DE genes but with the horizontal axis describing the amount of received signal. e Unique impact of CCC on downstream genes. A high score means the ligand-receptor pair is likely affecting the expression of the downstream gene considering the effects from other genes in the same cell.

Supplementary Figure 16
NOTCH signaling in human skin a Cell-cell communication directions. Sender: the direction shows to where the signal sending cells are sending the signal and color depicts the amount of signal sent by each cell. Receiver: the direction shows from where the signaling receiving cells are receiving the signal and color depicts the amount of signal received by each cell. b The cell-cell communication of all ligand-receptor pairs summarized to clusters (clustering identical to the recent NC paper), a thicker arrow means stronger signaling activity. c Dotplot showing the detailed signaling strength among the clusters through the individual ligand-receptor pairs. d The identified differentially expressed genes due to cell-cell communication activity. This is analogous to pseudotime DE genes but with the horizontal axis describing the amount of received signal. e Unique impact of CCC on downstream genes. A high score means the ligand-receptor pair is likely affecting the expression of the downstream gene considering the effects from other genes in the same cell.

Supplementary Figure 17
OXT CCC in MERFISH mouse hypothalamic preoptic region a Cell type plot of MERFISH mouse hypothalamic preoptic region data. b Heat map of cluster-level CCC where rows and columns correspond to senders and receivers, and the scatter plot of inferred sent and received signal through the ligand-receptor pair Oxt-Oxtr.

Supplementary Figure 18
CCC in STARmap mouse placenta data a Cell type plot of the STARmap mouse placenta data. b Cluster-level CCC heatmaps where rows and columns correspond to signal senders and receivers respectively and the spatial distribution of sent and received signal through the ligand-receptor pairs Igf2-Igfr2, Mdk-Sdc1, and Mdk-Sdc4.

Supplementary Figure 19
WNT signaling in seqFISH+ mouse secondary somatosensory cortex a Cell type plot of the seqFISH+ sscortex data. b Cluster-level CCC heatmaps where rows and columns correspond to signal senders and receivers respectively and the spatial distribution of sent and received signal through the ligand-receptor pairs from the WNT signaling pathway. PD1 signaling in the Visium data of human breast cancer.
a The amount of received signaling through two ligand-receptor pairs of PD1 signaling. b The signaling directions of the two ligand-receptor pairs.

Supplementary Figure 21
Cell-level correlation between inferred signaling and expression of known downstream genes Each dot represents a ligand-receptor pair. The Spearman's correlation coefficient is computed between the amount of received signal and the average expression of known downstream genes. The box plots show Q1 (25% quantile), Q2 (median), and Q3 (75% quantile). The upper whisker is the largest data point smaller than Q3+1.5IQR and the lower whisker is the smallest data point greater than Q1-1.5IQR. From left to right, n=192, 164, and 519 ligand-receptor pairs examined, respectively.

Supplementary Figure 22
Cluster-level correlation between inferred signaling and activity of known downstream genes (comparison with CellChat and Giotto) Each dot represents a ligand-receptor pair. The activity of known downstream genes of a ligandreceptor pair is quantified for each cell cluster as the percentage of positive significant differentially expressed ones (p-value < 0.05) in that cluster. Spearman's correlation coefficients are computed between the amount of received signal and the downstream gene activity. For COMMOT, the amount of received signal is quantified by the average received signal of the cells in that cluster. For Giotto and CellChat which output a cluster-level signaling table, the maximum CCC score with each cluster as the receiver is used to quantify the amount of received signal. The box plots show Q1 (25% quantile), Q2 (median), and Q3 (75% quantile). The upper whisker is the largest data point smaller than Q3+1.5IQR and the lower whisker is the smallest data point greater than Q1-1.5IQR. For COMMOT and CellChat results, n=119, 84, and 518 ligand-receptor pairs for the Visium human breast cancer, Visium mouse brain, and seqFISH+ datasets, respectively. For Giotto, the number of pairs are n=85, 79, and 339, respectively. The heteromeric ligand and receptors are converted to single unit pairs for Giotto and thus the number of pairs is different from the numbers in COMMOT and CellChat.

Supplementary Figure 23
Cluster-level correlation between inferred signaling and activity of known downstream genes (comparison with CellPhoneDB v3) Each dot represents a ligand-receptor pair. The activity of known downstream genes of a ligandreceptor pair is quantified for each cell cluster as the percentage of positive significant differentially expressed ones (p-value < 0.05) in that cluster. Spearman's correlation coefficients are computed between the amount of received signal and the downstream gene activity. For COMMOT, the amount of received signal is quantified by the average received signal of the cells in that cluster. For CellPhoneDB, the total CCC score with each cluster as the receiver is used to quantify the amount of received signal. For COMMOT, a uniform 1000 μm (spot/cell level) distance cutoff is applied. For CellPhoneDB, the distance between clusters is first quantified as the average distance between spots/cells of the two clusters, and different cutoffs are tried (5000 μm, 7500 μm, and 10000 μm for Visium datasets and 3000 μm, 5000 μm, and 7000 μm for the seqFISH+ data). The box plots show Q1 (25% quantile), Q2 (median), and Q3 (75% quantile). The upper whisker is the largest data point smaller than Q3+1.5IQR and the lower whisker is the smallest data point greater than Q1-1.5IQR. For the Visium human breast cancer, Visium mouse brain, and seqFISH+ datasets, n=119, 83, and 499 ligand-receptor pairs, respectively.

Supplementary Figure 24
Examples of inferred signaling and expression of known downstream genes.
The left column shows the ligands contributing to the inferred ligand-receptor complex. The middle column shows the inferred ligand-receptor complex. The right column shows the expression levels of known downstream target genes of the ligand-receptor pairs according to scSeqComm database.

Supplementary Figure 25
Comparison with CellChat on a Visium human breast cancer dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and CellChat, and by both methods. In comparison to CellChat, COMMOT uniquely detects CCC that are not differentially expressed in clusters but spatially significant.

Supplementary Figure 26
Comparison with CellChat on a Visium mouse brain dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and CellChat, and by both methods. In comparison to CellChat, COMMOT uniquely detects CCC that are not differentially expressed in clusters but spatially significant.

Supplementary Figure 27
Comparison with Giotto on a Visium human breast cancer dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and Giotto, and by both methods. Both COMMOT and Giotto prioritize CCC between spatially close clusters while COMMOT identifies more CCC links that are beyond the neighboring cells because Giotto uses cell-contact graph that might limit the range of inferred CCC. Giotto characterizes the average expression of ligands and receptors of spatially neighboring spots without providing the role for each spot in CCC. COMMOT identifies each spot as a sender, receiver or both in CCC with interpretable results.

Supplementary Figure 28
Comparison with Giotto on a Visium mouse brain dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and Giotto, and by both methods. Both COMMOT and Giotto prioritize CCC between spatially close clusters while COMMOT identifies more CCC links that are beyond the neighboring cells because Giotto uses cell-contact graph that might limit the range of inferred CCC. Giotto characterizes the average expression of ligands and receptors of spatially neighboring spots without providing the role for each spot in CCC. COMMOT identifies each spot as a sender, receiver or both in CCC with interpretable results.

Supplementary Figure 29
Comparison with CellPhoneDB v3 on a Visium human breast cancer dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and CellPhoneDB, and by both methods. CellPhoneDB v3 identifies more CCC links than COMMOT potentially due to its consideration of inter-cluster distances rather than distance among individual cells. COMMOT identifies CCC that is locally significant which CellPhoneDB v3 may neglect due to its consideration of expression level of entire clusters.

Supplementary Figure 30
Comparison with CellPhoneDB v3 on a Visium mouse brain dataset a The distribution of inter-cluster distances (distance between geometric center of the clusters) of the significant CCC cluster pairs. b Each point represent a ligand-receptor pair and a cluster pair. c-e Example significant cluster-level CCC that are uniquely identified by COMMOT and CellPhoneDB, and by both methods. CellPhoneDB v3 identifies more CCC links than COMMOT potentially due to its consideration of inter-cluster distances rather than distance among individual cells. COMMOT identifies CCC that is locally significant which CellPhoneDB v3 may neglect due to its consideration of expression level of entire clusters.

Supplementary Figure 31
Comparison of spot-level and cluster-level approaches on a Visium human breast cancer dataset a Clustering of the spatial data. b-c Left two columns: The level of sent and received signals of the spots inferred by COMMOT. Right: The significant ( 0.05) cluster-level signaling connections inferred by CellPhoneDB v3 with a cluster-level distance cutoff of 5000 µm. For the LR pairs ADIPOQ-ADIPOR1 and FGF7-FGFR1, COMMOT found local signaling hotspots (marked by dashed boxes) as subclusters of clusters 1 and 0, respectively, which were overlooked by the cluster-level inference.

Supplementary Figure 32
Comparison of spot-level and cluster-level approaches on a Visium mouse brain dataset a Clustering of the spatial data. b-c Left two columns: The level of sent and received signals of the spots inferred by COMMOT. Right: The significant ( 0.05) cluster-level signaling connections inferred by CellPhoneDB v3 with a cluster-level distance cutoff of 5000 µm. For the LR pair Adcyap1-Adcyapr1, COMMOT found a localized signaling hotspot within cluster 3 which is far away from other significant sender and receiver populations while the cluster-level method still identified significant signaling connections between cluster 3 and others. For the LR pair FGF7-FGFR1, COMMOT found local signaling hotspots which were overlooked by the cluster-level inference.

Supplementary Figure 33
Illustration of effect of competition on a Visium human breast cancer dataset a Clustering of the spatial data. b-c Top: The level of sent and received signals of the spots inferred by COMMOT with or without considering competition between ligand and receptor species. Bottom: The significant ( 0.05) cluster-level signaling connections inferred by COMMOT with or without considering competition and by CellPhoneDB v3 with a cluster-level distance cutoff of 5000 µm. For the LR pairs SPP1-ITGA5-ITGB1 and CXCL12-ACKR3, COMMOT with competition found localized signaling hotspots compared to the other two approaches.

Supplementary Figure 34
Illustration of effect of competition on a Visium mouse brain dataset a Clustering of the spatial data. b-c Top: The level of sent and received signals of the spots inferred by COMMOT with or without considering competition between ligand and receptor species. Bottom: The significant ( 0.05) cluster-level signaling connections inferred by COMMOT with or without considering competition and by CellPhoneDB v3 with a cluster-level distance cutoff of 5000 µm. For the LR pair Spp1-Itgav_Itgb1, COMMOT with considering competition found higher signaling activity on the right side of the tissue while the other approaches inferred active signaling everywhere. For the LR pair Igf1-Igf1r, COMMOT with considering competition also found localized signaling hotspots compared to other approaches.  1 Collective optimal transport in COMMOT

Collective optimal transport
Motivated by the biological application of analyzing interactions among multiple species of ligand and receptors, we introduce collective optimal transport (COT) which considers multiple source and target species, unlike traditional optimal transport [1,2] that only couples two single species. Consider n s points (cells of spots of cells) on which the discrete distributions are supported. Collective optimal transport couples n l source species (ligands in the COMMOT application), α i ∈ R ns + , i = 1, · · · , n l and n r target species (receptors in the COMMOT application), β j ∈ R ns + , j = 1, · · · , n r . Among these species, only part of the species pairs can interact and we use I to denote the indexes of iteracting species such that (i, j) ∈ I indicates that source species i interacts with target species j. For each species pair (i, j) ∈ I, a precomputed cost matrix C (i,j) ∈ R ns×ns + is given. An example of the cost matrix is the Euclidean distance between the n s points where a distance entry is replaced by ∞ if it exceeds the spatial range of signaling. Collective optimal transport finds a collection of coupling plans P (i,j) ∈ R ns×ns + , (i, j) ∈ I for each species pair in I that minimizes the total coupling cost over all pair.

1
The following optimization problem is formulated: where F penalizes the uncoupled masses µ i and ν j . In practice, there could be infinity entries in the cost matrices to block certain coupling according to specific interpretations in the applied domain. Also, α i and β j are not required to be probability distributions. This property enables the application of collective optimal transport to problems where the mass distribution should not be normalized, for example, when the units need to be preserved. Here, the units are the expression levels of various ligand and receptor species and they need not to be normalized individually so that they stay comparable.

Alternative formulation of COT
To solve the collective optimal transport problem described above, we reformulate the problem. The multiple source and target species are combined to formulate an optimal transport problem of size (n s * n l )×(n s * n r ). A new pair of distributions α ∈ R ns * n l and β ∈ R ns * nr to be coupled is generated such that α(i * n s + k) = α i (k) and β(j * n s + l) = β j (l) The corresponding new cost matrix is then defined as C(i * n l + 1 : (i + 1) * n l , j * n l + 1 : (j + 1) * n l ) = In this setup, the discrete OT problem is formulated as the following, where m = n s * n l , n = n s * n r , µ and ν are the unpaired masses, and F and G are penalty functions. In this work, we use two specific constructions with l 1 or l 2 penalty for the untransported mass. COT with l 1 penalty min P,µ,ν ⟨P, C⟩ F + ϵ p H(P) + ϵ µ H(µ) + ϵ ν H(ν) + ρ(∥µ∥ 1 + ∥ν∥ 1 ) s.t. P1 n = α − µ, COT with l 2 penalty min P,µ,ν In the equations above, H(x) = i x i (ln(x i )−1) is an entropy regularization function which is widely used to increase the computational efficiency [3]. The nonnegativity constraints in Eq. (3) are enforced here by the entropy barriers in Eqs. (4) and (5). The entropy terms also enforce a strong convexity which enables efficient numerical solutions.

Numerical solutions
In this section, we derive two kinds of numerical solutions to the optimization problems in Eqs. (4) and (5).
COT with l 1 penalty We begin with the Lagrangian associated with Eq. (4), The first order condition relates P with the multipliers as With (ϵH(·) + ρ∥ · ∥ 1 ) * (y) = ϵe y−ρ ϵ and condition in Eq. (7), we have the following Lagrange dual function, whose first order derivatives are These immediately leads to gradient based methods, such as momentum gradient ascent and Nesterov's method [4]. If we set ϵ = ϵ p = ϵ µ = ϵ ν , we can easily obtain the explicit first order conditions, We now apply the log-sum-exp trick to stabilize computation for small value of ϵ. Note that for any a ∈ R m . We now set a to f and g respectively and obtain a stabilized Sinkhorn iteration, for l ≥ 0 with arbitary f (0) and g (0) .
COT with l 2 penalty The algorithms for the COT with l 2 penalty can be derived similarly.
To avoid solving the highly nonlinear systems, we first introduce two primal variablesμ and ν leading to the following minimization problem, The associated Lagrangian then reads E(P, µ, ν,μ,ν; f , g, r, s) =⟨P, C⟩ F + ϵ p H(P) Taking the first order derivatives leads to the first order conditions, 4 We now obtain the Lagrange dual function by plugging these conditions to the Lagrangian E, L(f , g, r, s) = min P,µ,ν,μ,ν E(P, µ, ν,μ,ν; f , g, r, s) where (ϵH) * (y) = ϵe y ϵ and ( ρ 2 ∥ · ∥ 2 ) * (y) = ∥y∥ 2 2ρ . Similar to the l 1 case, gradient based optimization methods can be used given the first order derivatives, Setting ϵ = ϵ p = ϵ µ = ϵ ν , using the first order condition, and applying the log-sum-exp trick, we have the following Sinkhorn iteration, where ω is the Wright omega function.

Benchmark and comparison with PDE-based simulation
Ideally, CCC inference methods for spatial data should be validated by examining the spatial co-localization of ligand and receptor proteins. However, such data is often lacking. We therefore built partial differential equation (PDE) models to simulate CCC in space. The PDE model describes the essential physical processes in ligand-receptor interaction including the diffusion of ligands, binding and dissociation of ligand-receptor complexes, and degradation of molecules [5] in 1-dimensional (1D) or 2-dimensional (2D) spaces (Extended Data Fig. 1a).
As an example, in the case of two ligand species and one receptor species (Extended Data Fig. 1a, b), the two ligands could be Wnt2 and Wnt5a and the receptor could be Fzd5 [6]. Then [L 1 ] and [L 2 ] describe the spatial distribution of free Wnt2 and Wnt5a whose production rates are observed from spatial data. Similarly, [R] describes the spatial distribution of unoccupied Fzd5 whose initial distribution is taken from spatial data. Then [L 1 R] and [L 2 R] describe the spatial distribution of Wnt2-Fzd5 and Wnt5a-Fzd5 complexes, respectively, which reconstruct the CCC activity in space. The quality of the inferred CCC will be assessed based on how accurately it reproduces the simulated spatial distributions of the ligand-receptor complexes.

Model
Assume there are n l ligand and n r receptor species where each ligand or receptor species might bind to multiple others. Let [L i ], [R j ], [L i R j ] : R n grid ×2 → R denote the spatial distributions of ligand i, receptor j and the complex formed by their binding, the 1) diffusion and degradation of ligand and 2) binding and dissociation of ligand-receptor complex are modeled in the following equations, for (i, j) ∈ I, the index set of ligand and receptor that can bind. For each ligand-receptor pair, ligand i and receptor j that can bind, a (i,j) and b (i,j) are the association rate and dissociation rate between them. For ligand i, c i is its degradation rate and D i is its diffusivity in space.

Synthetic data generation
Synthetic benchmark data for validating the collective optimal transport method is generated by simulation of the PDE models described in Eq. (19). Various cases of ligandreceptor binding are considered which have varying numbers of ligand species, receptor species, and ligand-receptor complexes, to account for different complexities and levels of competition among species (Supplementary Figs. 1-4). For each case, ten independent simulations were carried out. For each simulation, the initial distribution of ligands and receptors were generated by using multiple Gaussian distributions with random bandwidths and random locations in a square domain. All ligand-receptor complexes were initiated to be zero.
For all simulations, a natural boundary condition (no-flux) was assumed and a 2D 50by-50 grid was used. The simulations were implemented using the py-pde package [7]. The solver was set to "scipy" which uses the "scipy.integrate.solve ivp" utility in SciPy [8] that implements the explicit Runge-Kutta method of order 5(4). The simulation was run until the concentrations of the ligand-receptor complexes reach (numerical) equilibrium. The concentrations of the ligand-receptor complexes depict the spatial distributions of each receptor species occupied by each ligand species. The data analysis methods (optimal transport based method) are then benchmarked based on their performance of reconstructing these spatial distributions.

Benchmark and comparison result
We first illustrate COMMOT in two simulated examples in 1D and 2D. Specifically, the inferred spatial distribution of bound ligand-receptor complexes is compared with the simulated one. In the 1D physical model, the two ligand species are located on the two sides of the location of receptor, respectively and will compete to bind to the receptor with a limit spatial range of interactions. The effects of symmetric ligand locations, biased ligand locations, and different ligand diffusivities are correctly captured by collective optimal transport (Extended Data Fig. 1a). In the 2D physical model, five interacting species are modeled, and the interaction patterns are also recapitulated (Extended Data Fig. 1b). In addition to revealing ligand-receptor complexes, collective optimal transport also annotates the source of ligands that contributed to the ligand-receptor complexes (Extended Data Fig. 1a,b).
Further quantitative comparison is carried out using simulations with randomly generated spatial expression of ligands and receptors in various cases of ligand-receptor binding (Extended Data Fig. 1c). To account to different ligand-receptor binding situations, ten different cases are considered with different numbers of ligand and receptor species and different levels of competition (the ratio of interacting pairs over total number of species) among the species (Extended Data Fig. 1d). For each case, ten random initializations are generated. The levels of bound ligand-receptor complex inferred by collective optimal transport are compared to the PDE simulation results using the root mean square error and Spearman's correlation coefficient. Our results show that collective optimal transport consistently derives accurate predictions across the different cases and outperforms the results by evaluating each ligand-receptor pair independently using traditional optimal transport. The improvement in the results from the pairwise approach is especially significant when competition among species is high (Extended Data Fig. 1d, Supplementary  Figs. 1-4), for example, cases 4, 5 and 10. A comparison with simulated data demonstrates that collective optimal transport used by COMMOT accurately reconstructs CCC in space. The ability to handle non-probability distribution is especially useful when the total levels of ligand and receptor expressions are significantly different (Supplementary Fig. 10). Also, incorporating spatial information helps prioritizing interactions among cells that are spatially close (Supplementary Fig. 11). Moreover, our method outperforms two other variants of optimal transport that also do not require the input data to be probability distributions: unbalanced optimal transport and partial optimal transport (Supplementary Figs. 5-6). In unbalanced optimal transport, the commonly used KL divergence is adapted in this comparison while other divergence forms may be incorporated as well. The differences among these optimal transport variants are further illustrated by both a one-dimensional example ( Supplementary Fig. 7) and with real ST data ( Supplementary Fig. 8-9), showing again the advantages in using COMMOT.

Joint analysis of mouse cortex with different datasets
We jointly analyzed CCC in mouse cortex datasets generated with three different technologies including Visium, seqFISH+, and STARmap (Extended Data Figs. [3][4][5]. Across these three datasets, AGT signaling was identified in neurons. Spatially, neurons in the L2-3 region were identified as strong receivers of AGT ligands across the three datasets. Interestingly, a striped signaling pattern was observed, wherein strong signals within individual layers form stripes, while weak signals form inter-stripe regions. Strong AGT signaling activity among oligodendrocytes was also identified in both STARmap and seqFISH+ datasets. (Extended Data Fig. 3).
In both Visium and seqFISH+ cortex datasets, we inferred WNT signaling to be active across different cortical layers. In both datasets, we identified WNT signaling to be relatively low in layer 5, compared to other layers. This result is in line with previous study showing that disruption to WNT signaling in Lrp6 mutant mice leads to cortical hypoplasia and reduced proliferation in neurons in layers 2-4 and 6, with neurons in layer 5 remaining relatively unchanged [9] (Extended Data Fig. 4).
TAC (tachykinin neuropeptide family) signaling activity was consistently found in both Visium and STARmap cortex datasets to be active in non-neuronal cells and in inhibitory neurons, especially in somatostatin-expressing neurons (Sst), because of the high expression of tachykinin receptor-3 in Sst neurons [10] (Extended Data Fig. 5).

Robustness evaluation of COMMOT on drosophila embryo data
For the robustness study, we utilized the stage 6 Drosophila embryo, an extensively studied system, which contains approximately 6000 cells that have been extensively studied with in situ databases available that present systematically annotated spatial gene expression [11,12]. The paired data consists of scRNA-seq of 1297 cells and 8925 genes along with the spatial single-cell resolution data of 3039 cells and 84 genes [13], with both good depth and spatial organization. We first integrated them using SpaOTsc [14] to generate two imputed datasets, the imputed scRNA-seq with predicted spatial coordinates of cells and the imputed spatial data with predicted expression of the 8925 genes from the scRNA-seq data. Using subsampling as a test for robustness, we randomly subsampled the data and compared them to the results from the full dataset. When decreasing the subsample size, the cosine distance (for spatial signaling direction) and Jaccard distance (for cluster level signaling) were found to remain relatively small showing that both signaling directions and cluster level signaling are robustly captured in the subsampled data (Extended Data Fig.  6). The differentially expressed signaling genes were then identified by using the inferred amount of received signal as the co-factor, analogous to the approach for pseudotime [15].

8
The DE genes identified from subsampled datasets were found to be comparable to the ones identified from the full dataset (Extended Data Fig. 6). The spatial signaling direction is described by a vector field defined on n grid points discretizing the tissue space, and is represented by an array V ∈ R n×2 . Cosine distance is used to compare the vector field V sub from subsampled data to the one from the full data V full and is defined as To compare two cluster level CCC networks S cl 1 and S cl 2 , we first binarize them such that the edges with p-value smaller than 0.05 are kept in the edge setsS cl 1 andS cl 2 . Then, the Jaccard distance is used for quantitative comparison, d Jaccard (S cl 1 ,S cl 2 ) = 1 − |S cl 1 ∩S cl 2 |/|S cl 1 ∪S cl 2 |.

Comparison to CellChat, Giotto, and CellphoneDB
Using Visium human breast cancer and mouse brain datasets, we studied the difference between COMMOT and the three other methods. When assessing the distributions of inter-cluster spatial distance of the identified significant CCCs, the significant CellChat interactions have a relative uniform distribution of the spatial distances. This is expected since CellChat was designed for scRNA-seq data and spatial organization is not considered in the method. Both COMMOT and Giotto prioritize CCC between spatially close clusters while COMMOT identifies more CCC links that are beyond the neighboring cells because Giotto uses cell-contact graph that might limit the range of inferred CCC. Cell-PhoneDB v3 identifies more CCC links than COMMOT potentially due to its consideration of inter-cluster distances rather than distance among individual cells . In comparison to CellChat, COMMOT uniquely detects CCC that are not differentially expressed in clusters but spatially significant (Supplementary Figs. 31-32). Giotto characterizes the average expression of ligands and receptors of spatially neighboring spots without providing the role for each spot in CCC. COMMOT identifies each spot as a sender, receiver or both in CCC with interpretable results (Supplementary Figs. 33-34). COMMOT identifies CCC that is locally significant which CellPhoneDB v3 may neglect due to its consideration of expression level of entire clusters ( Supplementary Figs. 35-36). Furthermore, COMMOT can identify localized signaling hotspots compared to clusterlevel approaches (Supplementary Figs. 37-38). For a specific ligand-receptor pair, COM-MOT further prioritizes regions containing its high signaling activity with low competitions from other pairs ( Supplementary Figs. 39-40).